


UNIVERS 
OF MICHIGAN 


JUL 2 1 1952 


ENGINEEA, 
LIBRARY ~ 








JULY, 1952 NuMBER 7 





CONTENTS 


On Weak Interaction of Strong Shock and Mach Waves Generated Downstream of 
Boa-Texs Cau 


Downwash Behind a Two-Dimensional Wing Oscillating in Plunging Motion 
E. LAPIN, R. CROOKSHANKS, AND H. F. HUNTER 
Analysis of Elastic Structures by Matrix Transformation with Special Regard to Semi- 
monocoque Structures BéryJE LANGEFORS 
Creep Buckling of Columns CHARLES LIBOVE 
Buckling of Sandwich Cylinders Under Bending and Combined Bending and Axial 
Cut-TEH WANG AND D. P. SULLIVAN 
The Effect of Nonuniform Surface Temperature on the Transient Aerodynamic Heat- 
ing of Thin-Skinned Bodies A. E. Bryson anp R. H. Epwarps 
Thermal Stress in Power-Producing Elements 
Visualization of Flow Fields by Use of a Tuft Grid Technique 
On Pohlhausen’s Method with Application to a Swirl Problem of Taylor. . J. C. Cooke 
Readers’ Forum: ‘The Stability of the Laminar Boundary Layer in a Compressible 
Fluid for the Case of Three-Dimensional Disturbances,” by D. W. Dunn and C. C. 
Lin; “Some Considerations of the Stability of Laminar Parallel Flows,’’ by Martin 
Lessen; “Taking Account of the Compressibility of the Material in the Plastic 
Buckling of Plates,” by P. P. Bijlaard, with ‘Authors’ Reply,” by Elbridge Z. Sto- 
well and Richard A. Pride; “‘A Generalized Ground Run Chart,” by William F. 
Savage; “Note on the Drag of a Finite Wedge at Mach Number 1,” by Isao 
Imai; ‘Note on the Conditions at a Sharp Leading Edge in a Supersonic Stream,” 
by Russell E. Duff; ‘Effect of Flapping Hinge Offset on Rotor Blade Stall,” by J. 
P. Chawla; “Comments on ‘Uniformly Loaded Semi-Infinite Wedge-Shaped 
Plates’—Author’s Reply,” by A. J. A. Morgan; “On the Boundary-Layer Equa- 
tions in Hypersonic Flow,” by S. F. Shen; “Pressure Distribution on an Airfoil in 
Nonuniform Motion,”’ by Robert H. Scanlan; ‘Analysis of the Elastic and Plastic 
Stability of Sandwich Plates by the Method of Split Rigidities—III,” by P. P. 
Bijlaard; ‘‘Once More—Single Degree of Freedom Flutter of an Aileron,” by K. P. 
Abichandani and R. M. Rosenberg 491-504 


Published by the 
INSTITUTE OF THE AERONAUTICAL SCIENCES, INC, 











JOURNAL OF THE 
AERONAUTICAL SCIENCES 





Editor: Hucu L. Drypen, Director, N.A.C.A. *- 
Associcte Editor: Ropert R. DExTER Managing Editor: Berneice H. Jaroge 





Editorial Committee 


Structures and Materials Flight Propulsion Instruments 


W. B. Bleakney George W. Brady 
E. W. Conlon Frank W. Caldwell 
L. H. Donnell L. S. Hobbs Charles H. Colvin 
E. C. Hartmann E. E. Johnson Dray 
N. J. Hoff Joseph Keenan 
R. P. Kroon 
Hans Reissner 
Abe Silverstein 
C. Fayette Taylor 
E. S. Taylor 
P. Taylor 
E. S. Thompson 
H. S. Tsien 


<rovn 


Chi-Teh Wang 
Aircraft Design 


Rotating Wing Aircraft 
Hall L. Hibbard 
Richard Hutton i 
Cc, L. Johnson R. L. ‘Bisplinghoff 
~ A. ee William Bollay 
Ralph MeClarren Frei E. Week 
RH. Prewitt Robert J. Woods 


Igor I. Sikorsky J. P. Den Hartog 


Physiologic Problems Fuels and Oils 
i D. P. Barnard 
J. H. Doolittle 
Otto E. Kirchner Graham Edgar 
Jerome Lederer R. T. Goodwin 
W. Randolph a II John C. Leslie S. D. Heron 
Ross A. MacFarland R. Dixon Speas A. M. Rothrock 





SUBSCRIPTION RATES 


OURNAL OF THE AERONAUTICAL SciENcES, Published Monthly.—United States and Possessions: 1 Year = 00; Single Copies,” i 

1.50. Foreign Countries Including Canada (American Currency Rates): 1 Year, $13.00; Single ‘Copies, $1.50 rah 
Fo rag Encineerinc Review, Published Monthly.—United States and Possessions: 1 Year, $3.00; Single Copies, $050, 

Foreign Countries Including Canada (American Currency Rates): 1 Year, $3.50; Single Copies, $0. 50. eb 
Notices of change of address should be sent to the Institute at least 30 days prior to actual change of address. 
Manuscripts for esteem proofs, and all correspondence should be addressed to the Editorial Office of the Journat. 
Correspondence regarding subscriptions may be addressed to Publication Office, 20th and Northampton Sts., Easton, Pa., or to a? 
Editorial Office of the Journa or THE AERONAUTICAL Sciences, 2 East 64th Street, New York 21, New York. 


Copyright, 1952, by the Institute of the Aeronautical Sciences, Inc. 


Entered as Second Class Matter at the Post Office, Easton, Pa., May 1, 1937. Acceptance for mailing at a special rate 2 
Postage provided for in the Act of August ‘24, 1912. Authorized April 29, 1937. te 














A 
shoc 
fron 


stuc 


dow 


flect 
duce 





JOURNAL OF THE 
AERONAUTICAL SCIENCES 


JULY, 1952 NUMBER 7 





VOLUM E 19 


—_—_ 


On Weak Interaction of Strong Shock and 
Mach Waves Generated Downstream 
of the Shock’ 


BOA-TEH CHU 
The Johns Hopkins Unwersity 





SUMMAR\ by a two-dimensional wedge in a supersonic stream. 
Now suppose that the profile of the wedge deviates 
shock wave produced by a wedge whose face is perturbed slightly | slightly from a straight line; there will be Mach waves 
from the straight-line configuration. The solution is used to generated at the surface of the wedge. The Mach 
study the interaction of a shock wave and Mach waves generated = wayes will interact with the shock so that the latter 
downstream of the shock. It is found that Mach waves are re ‘ : : , sala , ; 
also deviates slightly from its original straight-line 
flected from the shock wave without change of sign but with re . 
duced amplitude. The ‘reflection index” is given analytically 
It approaches zero as the shock wave degenerates into a Mach 


Analytical solutions are obtained for the flow behind an oblique 


configuration. As a result, the flow field behind the 
perturbed shock is no longer irrotational or isentropic. 
wave and increases monotonically as the shock strength to a Our problem is to determine the flow field behind the 
_ of ey ga than unity ae ae pbb shock wave. We are particularly interested in the na- 
shock wave alld the Mach waves generated Dy a sma yrotuber ° . - 

ance on the wedge surface siadies an “almost deta” tecia on ture of the reflection of Mach waves at the shock wave 
the shock wave with the longitudinal and transverse scales and the resulting change of shape of the latter, as well 
stretched in proper ratio. The shock wave also behaves like a 
porous elastic structure that deforms with a simple stress-strain the stream lines. 


as the entropy gradient and vorticity developed across 


telation. The solution is applied to calculate the pressure dis 
tribution on a “‘supersonic airfoil’’ at a large angle of attack. The 
: ;' ‘ é . . sea Method of Attack 
extent of the downstream propagation of a small disturbance 
eg., that caused by a small hump on the wedge surface—is also 
petimated wave is small, the fundamental equations governing 
the flow can be linearized. It is found that, in general, 
NTRODUCTION . ° : : 
INTRODUCTION the equations for the various flow variables assume the 


Since the disturbance produced behind the shock 


form of the simple wave equation with an entropy gra 

Statement of Problem é , I ys ‘ 
dient as a nonhomogeneous term in the equation (cf. 

I IS WELL KNOWN that, if the Mach Number of a reference 1), except those governing the perturbed pres 
supersonic stream is high enough, the stream, while — sure and flow direction which satisfy simply the homo 


flowing past a convex corner, will produce an attached geneous wave equation. As a result, it is possible to 
shock wave. The shock wave is straight and the flow write down the most general form of the pressure field 
behind it is uniform. For most cases the flow behind — and flow direction in terms of two arbitrary functions 
the shock wave is also supersonic. A simple applica say, F and Gand all others in terms of one more un 
tion of this flow configuration is the flow field produced = known function—the entropy distribution. However, 
one should not conclude that the pressure field and 


Presented at the Aerodynamics Session, Twentieth Annual 
Meeting, 1.A.S., January 28-31, February 1, 1952 
* The author wishes to take this opportunity to thank Dr. F 
Clauser for his helpful comments concerning this paper. To have to take into account the boundary conditions to 
Miss Vivian O’Brien, the author is indebted for her excellent i 
Brow; ; : be satisfied at the shock wave. 
rawings and assistance in preparing the final manuscript; to Ps oe _ 
Miss M. Ann Emmart, for her help in typing up the manuscript. rhe boundary conditions to be satisfied are the tan 
t Graduate Student, Department of Aeronautics 


flow direction are therefore independent of the rota 
tionality and anisentropic nature of the field, since we 


gency of flow at the surface of the wedge and the con 
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tinuity, momentum, and energy conditions across the 
shock wave. The energy condition is automatically 
satisfied by the perturbed field, since the total tem- 
perature is constant everywhere in the field~-ahead of 


the shock as well as behind it. The condition of tan- 
gency at the surface of the wedge and one continuity 
relation and two the 
shock provide four boundary conditions from which the 


functions (/, G), the entropy distribu- 


momentum conditions across 
two arbitrary 
tion, and the shock angle can be determined. These 
four equations can be combined into a functional equa- 
tion of F. 
responding to the present problem is given. 
eral solution of the functional equation is also ob- 
The flow field behind the shock wave is thus 


The solution of the functional equation cor- 
The gen- 


tained. 
completely given in analytical form. 


Previous Related Work* 


Step-by-step calculation of the interaction problem 
of a shock wave and Mach waves by means of the char- 
acteristics method is well known (e.g., references 2 and 
3). For analytical treatment, simplifications are ob- 
tained by restricting ourselves to some limiting cases. 
The first type of simplification is based on ‘‘weak inter- 
action’’—that is, the shock wave is only slightly per- 
turbed from its normal configuration as a result of the 
interaction with Mach waves. Two of the previous 
investigations belonging to this category are those of 
Burgers‘ and Adams.° In fact, the author is greatly 
indebted to Burgers’ excellent paper, which inspired 
the present investigation. Burgers studied the trans- 
mission of sound waves through a shock wave in a one- 
dimensional flow. Adams, using the differential equa- 
tions developed by Sears! investigated the interaction 
of a normal shock wave and Mach waves generated 
ahead of the shock. A similar idea has been employed 
by Lighthill® * and Ting and Ludloff* in their study of 
the diffraction of shock waves. In all these analyses, 
the shock wave is not necessarily weak; therefore en- 
tropy change and rotationality downstream of the 
shock wave have to be accounted for. Since the shock 
wave is only slightly perturbed from its normal con- 
figuration, linearization of the differential equations 
governing the flow downstream of the shock wave is 
As the title indicates, the present paper also 
In the title, the adjective 


possible. 
belongs in this category. 
“‘strong’’ has been used to describe the ‘‘shock’’ for em- 
phasizing that the latter need not be weak. For 
“tweak shocks,”’ another better simplification is avail- 


able. This type of simplification is obtained as a re- 


* After this paper had been presented, my attention was called 
to a similar work by S. I. Pai, Cornell Aeronautical Laboratory, 


Inc. Through communication with Dr. Pai, it was learned that 
he had worked out the present case, as well as the axially sym 
metric case with a different method, two years ago. Unfortu 
nately, however, Dr. Pai’s report is still unavailable; hence, a 
comparison of the two studies cannot be made. Dr. Pai assured 


the writer that the paper will be made available in the future 
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sult of neglecting terms of higher than second order in 


the shock strength. 


As a consequence, there wil] be 


no change in entropy across the shock wave and no re. 


flection of waves at the shock, and the characteristic 


theory of the irrotational flow can be applied. 


This 


idea was employed by Chandrasekhar® (and later py 


others) to study a special interaction problem. The 


complete theory of weak shock was developed py 


Courant and Friedrichs. " 


It is interesting to note that 


Busemann'' has probably applied the idea some 2 


years ago in determining the shape of the shock wave 


produced by a flat plate at an angle of attack. 


Fried 


richs'* applied the method to study the formation and 


decay of shock, as well as the asymptotic nature of the 


shock 


wave 


in two-dimensional steady flow. The 


two-dimensional case was discussed earlier by Light. 


b,*? 


was investigated by Whitham.'' 


The analogous problem for a body of revolution 


In all these anal 


yses, aside from assuming that the shock wave is 


weak or only moderate so that entropy change and 


reflection at the shock may be neglected, no further 


simplifying assumption is made with regard to the flow 


downstream of the shock. 


The problem of the forma- 


tion and growth of a shock wave was independently 


tackled by Pillow,'® who also gave a solution account 


ing for the third-order terms in the shock strength. 


Aly) 
Ge i 


f(x) 
F(x), 
M, 


G(x) 


NOTATION 


an arbitrary function of y 

specific heats at constant pressure and volume, 
respectively 

equation of the profile, f’(v) df(x)/dx 

arbitrary functions 

Mach 


oblique shock 


Number of flow behind the undisturbed 
pressure 

gas constant 

entropy 

temperature 

stagnation temperature 
velocity components in the x- and y-direction 
Cartesian coordinates 

angle between the x-axis and oncoming wave 
When the later taken 
along the wedge surface, a is the half angl 


(Fig. 1) Y-axis is 


of the wedge 


V Mi? I 

Gost 

sin (a; + @)/sin (yy w)] ( l 

defined by Eq. (44 

density 

the “reflection index,” |sin (A wy) /sin (A 
Mi) | \ l 


Mach angle 

w Q 

wave angle 

w when the x-axis is taken along the wedge It 


is therefore the angle between the unper 
turbed shock wave and the direction of flow 
immediately behind it 


vorticity 
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Subscripts 
= condition of flow ahead of the shock, except Xp, 


which 


represents a point on the wedge, and Fi(x x which re pre 
sents a particular solution of the functional equation 


Eq. (A-1) | 
= condition of flow behind the unperturbed shock 


Tur DIFFERENTIAL EQUATIONS AND THEIR GENERAL 
SOLUTIONS 


Ahead of the shock wave, the stream is uniform and 
undisturbed. Downstream of the shock wave, the flow 


is governed by the following set of equations: 


Continuity Equation 


(Opu/Ox) + (Op7/Oyv) = 0 (1) 
Momentum Equations 
Ou Ou | Op 
u — v = — (2) 
Ox oy p OX 
Ov Ov 1 Op 
u rv =— (33) 
Ov Oy p Oy 
Energy Equation 
Cel + (1/2)u" + (1/2)2? = Cote (4) 
Gas Law 
p - pRT (5) 


The five equations are sufficient to determine the five 
unknowns p, p, 7’, “, and v—the pressure, density, tem- 
perature, and the x- and y-components of the velocity. 
It is convenient to choose the x-axis along the wedge 
surface and the y-axis perpendicular to it. 

A direct consequence of the energy and the momen 
tum equations is that downstream of the shock the 
entropy is constant along each stream line. For, by 


eliminating w*, v? from Eqs. (2), (3), and (4), we have 


u (O0S/Ox) + v(O0S/dOyv) = 0 (6) 


where the entropy S is given by 


S— S C, log (p/ pi) (pi/p)’ (7) 
The subscript | in the above equation represents some 
reference state that will be chosen later. 

Eq. (6) can be used to replace one of the three equa 
tions |Eqs. (2), (3), and (4)| from which it is derived. 
We shall replace Eq. (2) by Eqs. (6) and (7). 

If the wedge profile is strictly a straight line, then the 
flow behind the shock wave is uniform and irrotational. 
The shock wave is straight, and the entropy is the same 
for all points downstream of the shock wave. Let us 
denote the pressure, density, temperature, entropy, and 
velocity components parallel and perpendicular to the 
wedge surface behind the shock wave by /, pi, 71, Si, 
m1, %) (= 0). If the wedge profile deviates from a 
Straight-line profile, the disturbance generated at the 
boundary will interact with the shock wave, so that 
the latter also deviates from its normal configuration 
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and becomes curved. As a result, the flow behind the 
shock wave is no longer irrotational or isentropic. 
However, when the deviation of wedge profile from its 
normal straight-line profile is small, the shock wave will 
also deviate only slightly from its normal configuration, 
and the variation of p, p, 7, S, u, and v from their nor- 
mal values /;, pi, 71, #1, and v7, (= 0) will also be small. 
Consequently, one may put 

p = pit op 

p = pi t+ dp 
etc., where 6p/ pi, 6p/pi, 67°/7), S/R, bu/m, 60/u, < 1. 
Substituting these into Eqs. (1), (3), (4), (5), (6), and 
(7) and retaining only small quantities of first order, 


O (°*) Oo (~) Oo (*) 
+ + = () (S) 
Ov \pi Ov \uy Ov \m 
O (“") l oO (°") 9 
_ (9) 
Ox \u yA1\? Oy \pi 
l 6T° bu 
(= + us ( = 0 (10) 
| fe l Ti) iy 


6p Pr = (dp pi) + (67 7) (11) 
(0/Ox) (6S/R) = 0 (12) 
(6.8/C,) (6p/ pi) — y(6p/pr (13) 


Thus we obtain a simple linear system of six equations 
5p/ pi, 6p/p1, 67/7), 6S R, 6u/u4, 
namely, Eqs. (10), 
(11), and (13) More- 


over, Eq. (12) can be integrated immediately, giving 


and six unknowns: 
6v/u;. Three of the six equations 


are already in algebraic form. 


6S/R = A{y) (12’) 
that is, 6S is a function of v only. 

To see the effect of entropy gradient on the various 
flow parameters, we shall derive a single differential 


equation for each flow parameter. Thus, by elimi- 


. o? (sp Oo? (bp 
(Mi? — 1) — - ( =0 (14) 
Ox* \pi Oy \pi 
O° /év O° / dv 
(M,? — 1) ( ) - ( ) = 0 (15) 
Ox? (My Oy> \m 
af O° /6 fr s85 
(me — 1) 2 (°) : (”) (=) (16) 
Ox" Pi Oy Pi dy’ Us 


nation, 


oO” /éu O° /éu 
(\7,7 -— 1) 4 —= é 
Ox? \iy Ov- \m) 
P ~) 
( (1S) 
yMidy? \R 
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First of all, we see that, when /, > 1, the differential 
equations are of the type of a simple wave equation.* 
It is interesting that the equations governing the per- 
turbed density, temperature, and the x-components of 
the perturbed velocity involve a nonhomogeneous term 
containing the entropy change. On the other hand, the 
differential equations governing 6p/p; and 6v/u are 
not affected by the fact that the flow is slightly rota- 
tional. However, this does not mean that the pressure 
field and flow direction are independent of the rota- 
tionality and anisentropy of the field, because we have 
to take into account the boundary conditions at the 
shock wave. These boundary conditions will be dis- 
cussed later. 

Aside from illustrating the effect of entropy distribu- 
tion on the form of the differential equation, thus mak- 
ing possible a wise choice of the differential equation to 
be used in any particular problem of this type, there is 
actually no necessity of writing out, as we have done, 
all of the five equations {[Eqs. (14)—(18)], for the solu- 
tion of any one of them determines completely the 
values of the other unknowns by the set of relations 
(8), (9), (10), (11), (12), and (13). For our problem we 


shall use Eq. (15). The general solution of it is well 


known, 

év/u, = F(x — Biy) + G(x + Bry) (20) 
where F and G are two arbitrary functions, 6, = 
VM — 1. Then it follows from Eqs. (8)—(13) that 

Spr 

pyM’? 

; [F(x — By) ~ G(x + Bry) | (21) 

VM, - 1 

dp M,? 


= [F(x — Bw) — 
Pi V M;? — | 


5S 
G(x + Biy)] — ~ (22) 


p 


ST (y¥- DM? 


a [Fix — fay) — 
T Vv M,? — | 


G(x + Biy)| + (23) 


p 


*When 1, these equations become either Laplace or 


Poisson equations. 
duced to a system of two equations and two unknowns which may 


The system of six equations can also be re- 


be particularly useful in case 7, <1. Thus, we have 


re) (*") ; l M,2 Oo (°?) 
Oy \m 7 yMi? ox \ pi 
re) bv l re) 6 
("= ag(2) am 
Ox \ uy yM)\? oy \pi 


which may be reduced to the Cauchy-Riemann equations with 
the dependent variables 6v/m, 6pV 1 M,?/piyM\? and inde 


pendent variables x/V 1 M\?, yor x, yV 1 — M,?. 
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bu l 
My Vv 





Fix — By) — 
G(x + Biy) | = : (24 
2p & 


The vorticity distribution can be easily calculated, for 
by definition, 


Ov Ou O o .. 
c= = = (dv) — (du) (95 
Ox oy Ox Oy 
forourcase. Hence, 
¢ = (u/yMM,") (d/fy) (6S/R) (26 


which could have been derived directly from Eqs. (9 
and (10) or from the Bjerknes formula.t 
that the vorticity is also a function of y alone. 


It is seen 


BOUNDARY CONDITIONS 


The first boundary condition that has to be satisfied is 
that the flow must be tangential to the wedge surface 
y=f(x). Thatis, 


dv/u, = dy/dx = f’(x) (27 


where the relation is to be applied at y = 0. The other 
boundary conditions come from the shock wave. In 
short, the conservation of mass, momentumf and en 
ergy of the fluid crossing the shock wave must be en 
sured. 
wave by subscript 0, these conditions are (Fig. 1) 


Denoting the condition ahead of the shock 


polo SIN wo = p(u sin w — Vv COS w) (28 
Uy COS W = UCOSw + Vv sin w (29 

Po + pylty” sin? > 
P + potty SiN wy (U Sinw — vecosw) (30 


Cero = (1 2) uo" = Cot + (1 2)u? + (1 2)y" (31 
Po polo = p pl —= R (32 


where w = a» a is introduced for brevity, a being 
the angle between the oncoming stream ahead of the 
shock wave and the x-axis. 

It is well known that, if the wedge has a strictly 
straight-line profile, the shock 
straight and that the flow behind it is uniform and 1s 


wave will also be 


obtained from Eqs. (28)—(32) by putting vy = 0. That 
iS, 

polly SIN wo = pil; SIN w (33 

Uy COS Wy = Uy COS a (34 

Po + potty” sin? wo = Py + potty SiN wy(M, Sin a) (35 


Sears.' By a perturbed 


stream function and making use of Eq. (26), it can be readily ver 


+ This was done by introducing 
fied that Eq. (15) can be converted to the equation derived by 
Sears. 

t Used in the sense that pressure is considered as an equivalent 


form of momentum. 





whi 
diffe 
u,v, 
tion 
54) 
5p, 
as t 
conc 
field 
It 
4]) 
with 
beca 
Equi 
the ( 
satis 
grad 
thert 
surfa 
tion 
the « 


(24) 


d, for, 


(26 
{S. (9) 


S seen 


fied is 
urface 


(27 


other 

In 
d en- 
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(29 
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CyTo + (1/2)t0? = Cpt + (1/2)u,? (36) 
Po poly = pi/mli = R (37) 


from which pi, pi, 71, %1, and w can be calculated (as 


— a@and a = constant here). These are, in 


wy = WO 


fact, the usual ‘shock relations.” 

Now, if the wedge profile deviates slightly from a 
straight-line profile, the shock wave will also curve 
slightly about its normal position y = x tan a, (see 
Fig. 1), and the wave angle w) will vary slightly along 
the shock wave. At any particular point on the shock 
wave, let the change in the wave angle «» be dw; then, 
corresponding to this slight change of the wave angle, 
the condition of the flow p, p, 7, u, v behind the shock 
wave will also change slightly to p, + 6p, pi: + dp, 
T, + 67, u; + du, + 6v. The increments 6), dp, etc., 


are related by 


6p : bu 2 6v 
sin w) + sin w, — COS w, + 
pi iy 1 
Po a 
(: — ) cos w1:dw = O (38) 
Pi 
bu év\ ; 
COS @)1 T sin @, + 
Uy i) 
Pl Po n 2 
(: — ) SIN w,: day = O 
Po Pi 


bp éu\ . , év\ 
~ sin? w, — SiN w; COS wr — 
bryMy? My Uy 


Po . 
(1 _ ) SIN @; COS w:dw@ = 0 (40) 


39) 


a 


P1 
l 6T Ou 
(5. +( M;* = 0 (41) 
7 l T; uy 
6p/ pi = (6p/p1) + (67/7T)) (42) 


which may be derived from Eqs. (28)—(32) by taking 
differentials of these equations, considering p, p, 7, 
u,v, and w as tunctions of w). In deriving these equa- 
tions, use has been made of the equalities | Eqs. (33) and 
(34) | so as to simplify the resulting expressions. Since 
5p, dp, etc., introduced here have the same meaning 
as those introduced in the last section, these are the 
conditions that must be satisfied by the perturbed 
field at the shock wave y = x tan a. 

It is obvious that the perturbed field satisfies Eqs. 
41) and (42) automatically because they are identical 
with Eqs. (10) and (11). That this is so is no accident, 
because these equations are derived from the Energy 
Equation [C,7 + (1/2)u2+ (1/2)v? = constant] and 
the Gas Law (p = pRT), respectively, both of which are 
satisfied everywhere in the field. Should there be any 
gradient in the stagnation temperature, either due to a 
thermal discontinuity —e.g., a flame front or a contact 
surface—or to an initial nonuniformity in the stagna- 
tion temperature—e.g., the aniso-energic shear flow 
the energy equation [Eq. (41)] will form an additional 
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w, = wave angle 











direction of oncoming flow 


Fic. 1. Derivation of boundary conditions at the shock wave 
Equation of the unperturbed shock wave: y = x tan aw;. Dotted 
line indicates portion of the perturbed shock. The perturbed 
velocity components 6u, dv have not been shown in the figure 


equation from which the distribution of the stagnation 
temperature can be determined. 

Returning now to the problem of the flow behind a 
slightly perturbed shock wave, the three conditions 
that must be satisfied by the perturbed flow are given 
by Eqs. (38), (39), and (40). We may call them the 
conditions at the shock Conse- 


boundary wave. 


quently, we have now four boundary conditions 
three at the shock and one at the surface of the wedge 
to be satisfied. We have shown in last section that 
the perturbed flow field can be given in terms of three 
arbitrary functions—-namely, F(x — fy), G(x + 
By) and 6S(y). If we substitute Eqs. (20), (21), (22), 
(23), and (24) into the four boundary conditions | Eqs. 
(27), (38), (89), and (40) ], there will be four unknowns 
F, G, 6S, and the local wave angle dw). Thus the four 
boundary conditions are just enough for us to deter- 
mine the four unknowns. 

Instead of solving simultaneously the four equations 
for the four unknowns, we shall attempt to simplify 
the analysis by reducing the system to one of solving 
Such 
a reduction is desirable, especially because we have 
seen that the perturbed pressure field 6p/piy./,"° and 
the flow direction 6v/u; do not contain 6S explicitly and 
Now the boundary 
It is clear, 


a system of two equations and two unknowns. 


involve only the functions F and G. 
condition [Eq. (27)] involves only 60/1. 
therefore, that, if we can reduce the other three bound- 
ary conditions into a relation between 6p/piy./,? and 
6v/u,, then the two functions F and G can be deter- 
mined immediately. Once F and G are found, the en- 
tropy 6S and the local wave angle da» can be calculated 
from the remaining two boundary conditions at the 
shock wave. 

To obtain a relation between 6p/piyM,? and 6v/u, at 
the shock wave, we need only eliminate 6p/p,, 67/74, 
6u/%, bw) from the five equations, Eqs. (38)—(42). 
The result of this elimination can be conveniently 
written as 
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6p/piyM? = (6v/u) tan A 12 


where* 


ee tet 
( we 
: ) — , * sin? > sin” y - . > 
(1 ~ Me l + - M7? sin: w1 } + {| y.1/,? sin? a, — a vo 


tan A = 


* The formula for tan \ can be written in various forms. 


tan \ = sin? wy sin 2w;/(sin? 


where yo, wi are the Mach angles ahead and behind the shock wave 


is a function of J/; and w,, the values of which are 
known. In the derivation of Eqs. (48) and (44), use 


has been made of the shock-wave relation 


es... 2 | Y | 
po y + 1 A? sin? @ y+ 1 


The boundary condition [Eq. (43)]| is to be applied at 


y = x tan w —1.e., the undisturbed position of the shock 
wave. 

Summarizing, the two boundary conditions from 
which F(x — @1y), G(x + B1y) can be evaluated are: 
aty = 0, 


60/%, = f(x) (27 


at y = x tan a), 


6p/priy Mi? = (6v/im1) tan d» (43 


Once the functions F and G are determined, we can 
calculate the wave angle 6w) and the entropy 6S from 
39), and 


any two of the three equations, Eqs. (38), 


(40). In fact, it is simple to show that, at y = x tan 
Wy 
, : oe. I ép = 
dw 2) 


2 sin 20, pryM;" 
6S x ( l = M,? sin” W@W) ‘+ x 
R 2, ™ . 4 er 
sin? a (1 + 2 M,? sin? w, 


6p Pry M,? (46) 


Since it is essential in the following discussion to 
know the magnitude of \ as given by Eq. (44), we shall 
devote a paragraph to examining this equation. We 
notice that, when the shock wave degenerates into a 


Mach wave, 


tan A > tan w, = 1/* M,? — | 
so that 
6p l l év 
pyMP VMP2-—1% 


which is the usual characteristic relation. On the other 


hand, tan \ = © when the denominator 


For = 


wy COS” yb 


A\/7,? sin? «1) sin 2w; 


xample, 


S1N* wo COS* a ) 


l “Poe, ae 
(1 we : (1 + M; SIN” @, 
M,, 2 
(van sin? w . a - ) cos 


vanishes. This may be shown to be the correct rela 
tion connecting .J/; and w; when the shock wave assumes 
the position of maximum flow deflection. In fact, 
this can also be seen directly from Eq. (43); for, when 
tanA = © orA = 7/2, bv/m% 0 for a finite increment 
of pressure—a condition that can occur only when the 
shock wave is at the position of maximum flow deflec 
tion. 

The variation of tan \ with J/; at constant «@ is 
sketched in Fig. 2. (Variation of a, will only slightly 


alter the position of the curve.) The ordinate assumes 


the value tan yu; (the tangent of the Mach angle—i.e,, 
1/V M, 1) at the right end of the curve at which the 
shock wave degenerates into a Mach wave. As J 
decreases, tan A increases monotonically. It shoots 


up to positive infinity when /; is slightly less than 
unity, at which the shock wave assumes the position 
of maximum flow deflection. It becomes negativel 
infinite when .J/; is further reduced. However, for 
our problem, this part of the curve is not used and, 
hence, is not shown in the figure. In our case, there 


fore, 


This inequality will be useful in later analyses. 
When Eqs. (20) and (21) are substituted into Eqs 


27) and (43), we obtain a pair of simultaneous equa 
tions for FandG. If we eliminate G from the two equa 
tions, we obtain a functional equation for F. The 
general solution of the functional equation contaims 
another arbitrary function that has to be evaluated 
from a fifth boundary condition. This fifth boundary 
condition depends upon the type of flow behind the 
shock wave. For our problem, if the flow behind the 
shock wave is subsonic, the fifth boundary condition 
appears as the requirement that the disturbances must 
die down at infinity. On the other hand, if the flow 
behind the shock wave is supersonic, the fifth boundary 


condition appears as certain regularity requirements 
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that have to be satisfied at the “tip” of the body. 
More explicitly, the solution has to remain finite as x 
spproaches zero (from the right). 

REFLECTION AT THE SHOCK WAVE 


The two equations governing the F and G functions 


are: atv Q, 
6v/u, = f’(x) (27) 

aty = V tana, 
bp pry? = (6v/m1) tan (43) 


Now we have shown that the general solutions of the 
differential equations describing the flow field behind 


the shock wave are given by Eqs. (20)—(24). Substi- 
tute Eqs. (20) and (21) in Eqs. (27) and (43), 
F(x) + G(x) = f’(x) (47) 
tan A tan w)Fl[x(1 B, tan w,)] + 
tan A + tan w,)G[x(1 + 6, tan w,)] = 0 (48) 
where tan wn; = | V VW)? — 1 = 1/8, or wy = the Mach 
angle. The first equation clearly governs the reflec- 


tion at the solid surface, while the second one governs 
the reflection at the shock wave. 

Eliminating the function G from these equations, 
we obtain a functional equation for F, 


tan AX — tan yp, 


F Y l T Dy 
tan A + tan yp; 


tan w)) 


Fi[x(1 — 6; tan w,)] = f’[x(1 + B, tan a) (49a 


which mav also be written as 


L tan yw; + tan a 
|: 
tan uw; — tan wa, 


Fix i 


tan A — tan yp; 


tan A + tan yu; 
tan yw; + tan a, 

v (49b) 
tan w 


| tan py 


“ ) sin (A uy x 


sin (A + py 


F(x fr (« SIn (41 +7 ) 19¢ 


sin (iy Wy 


, Sin (pu 
I (» 
Sin (74 


Che last form is the most convenient one for the ensuing 
discussions. 

Some information concerning the nature of the re 
flection at the shock wave and the flow field behind it 
can be obtained without solving the functional equa 
tion Eq. (49c). The complete solution of Eq. (49 
:pplied to our problem will be given in next section. 

In Eq. (49¢), /’(x) is a known function of x. The 
equation shows that, if we are given the value of F(x 
ata certain point x, then we can calculate the value of 
, where 


F(x) at the points xT’, xT, xT, ... 
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tan 






«),= Constant 





tan = tan«é =——— 
2 
\/ MF -1 
M, =| M, =——— | 
1 sin @, 

Fic. 2. Variation of tan \ with 1/, at constant a. The posi 
tion of the curve is slightly changed for different values of w 
When the shock wave degenerates into a Mach wave, \ = wy 
When the shock wave is at the position of maximum flow deflec 
tion, tanA = « For MM, > 1,7/2 >A> wm > w. 
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S f yp 
P 
eT 
ota X(=Xo) 
Fic. 3. The wave produced at xo hits the wedge again at xol 
xT? 
Cr = sin (uy + w)*/sin (yy Ww} (50 
This suggests that the points xT, xl’, xT’, , may 
be the first, second, third, ... point of reflection on the 


solid surface of a disturbance generated at x as shown 
in Fig. 3. 

To verify the above speculation, we shall assume 
that a small disturbance is introduced at the point x 
This disturbance, according to the 
Eqs. (14)—(18 
tan yu, until it hits 


Xo (see Fig. 3). 
differential equations , will propagate 
along the Mach line y x Xo 
at which it is reflected 
Mach The 
point of intersection of the Mach wave and the shock 


the shock wave y 
back along the other family of 


x tan w, 


waves. 


wave 1s clearly 


\ Yo |tan w,/ (tan py tan w, dla 


y Xp [tan vw, tan w)/(tan py tan w Slb 


The reflected wave will therefore have the equation 


* Note that uw; > for a shock wave of finite strength so that, 
in general, [ l I’ can also be written as 
I (tan gw, + tan w tan wu tan 
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Thus, indeed, the disturbance generated at the point ie 
xo is reflected back by the shock wave, hits the wedge - 

° 7 " z 5 de 
surface at the point xI, and is then reflected away poi 

y from the surface and then back again from the Shock Eq 
Ss wave, hitting the wedge at the point xI*, ete, p hoc 

. S. 

k (> 1), therefore, has the significance of being the ratio Con: 
of the distances measured from the leading edge of ty ene 
successwe points of reflection on the wedge surface for « ate 

7 re a a 
given disturbance. From Eq. (50), it is clear that the dam; 
value of IT’ will increase if the shock strength jis de. also 
“LLL / creased, and in the limit when the shock wave degen. show 
Fic. 4. Starting from a point (x, x tan w,) on the shock, the crates into a Mach pana [ —+ o. When the flow TI 
forward propagating Mach wave and backward propagating behind the shock wave is subsonic, I’ becomes imagi be W! 
Mach wave intersect the wedge surface at x(1 — §; tan w;) and nary so that no reflection phenomena will be obe : 
x(1 + B, tan w), respectively. my ‘ - 10n phenomena will be observed 
One may also easily find the physical significance oj 
the quantities x(1 — 6; tan w) and x(1 + 38; tang from 
involved in Eq. (49a), as well as in subsequent discus. wave 
sions. It is clear from Eq. (5la) that lectit 
Ww 
tan uw; — tan a, flo 
Xx =x - = x(1 — 2, tan a) shock 
an py , 
flecti 
That is, if the abscissa of a point on the shock wave is no lo 
x, then the Mach wave passing through this point and Fit 
propagating forward will hit the surface at a distance wave 
x(1 — £; tan w) from the tip of the wedge. Similarly, cordi 
the Mach wave passing through the point on the shock 
wave and propagating downstream will hit the surface 
at x(1 + 6, tan a). These are shown in Fig. 4. 
SL / Next, suppose that the initial disturbance as given A ph 
Fic. 5. Compression wave is reflected from the shock wave asa _ by /"(x) 1s limited to a finite distance from the tip of the one 
weaker compression wave. wedge, say, x <x). Then, for x > 1, f(x) = 0 and the pressi 
functional equation [Eq. (49c)| reduces to, for x > x, the s 
share : : then 
F(xT) — oF(x) = O (52 
stress 
where 2/(¥ 
’ ; - dictir 
o = sin (A — pwy)/sin (A + wy (53 
imme 
Note that 0 < o < 1, because 7/2 > X > yw; for our case versa 
It follows immediately from Eq. (52) that, for x > x, hind 
eee : for, a 
F(xl) = oF (x) bely 
yee nee ehin 
FP(x¥?) = oF (x) , 
(94 
F(x¥") = o F(x) pry: 
~ which demonstrates that the disturbance becomes 
Fic. 6. The boundary casts a similar image on the shock weaker and weaker as we go downstream (since ¢ < | - 
wave except that the longitudinal and transverse scale are : ; kt Using 
stretched by some factor. The quantity o as defined by Eq. (53) behaves as the rt- 
wave, 
flection index or as the damping factor per cycle of tt 
‘tio Since it is well k y ‘ flecti F Mach 5p 
tan py tan a flection. Since it is well known that reflection of M 6f 
y— *ee oles tan (t — wi) X waves at solid surface does not alter the strength and pry] 
an uw; — tan a, : : . : 
sign of the wave, the damping must therefore occur at 
tan ‘ ; ‘ F 
E =" nm the shock wave. In fact, this may be seen directly 
tan wu; — tan w from Eq. (48). Thus, solving for G and introducing so tha 
“m™.° : : , . V. as de > r Ras. (5 ‘ 53 ) &) be- 
This reflection wave will hit the surface of the wedge ! and o as defined by Egs. (50) and (53), Eq. (4 ; 
i comes 00) a 
y = Oat the point 
G(xl) = —oF(x) (59 
tan wi + tan o r 7 
x = Xo = X% : . ~~ 
tan uw, — tan a, which states that the incoming wave at the point «1 lhese 
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These relations are valid in general. 
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ie. the wave once reflected from the shock wave) is 
damped down from the incident wave generated at the 
point x by the factor ¢. By substituting Eq. (55) into 
Eq. (21), it is clear that the reflected wave from the 
shock wave has the same sign as the incident wave. 
Consequently, a compression (or expansion) Mach wave 
venerated downstream of a shock is reflected at the shock 
wane as another compression (or expansion) wave with 
damped amplitude. Without using Eq. (21), this can 
also be readily seen from the geometrical construction 
shown in Fig. 5. 

The reflection index as defined by Eq. (53) can also 
be written in the form 

o = (tanA — tan w),)/(tanA + tan pw) (53a) 

from which it becomes evident that the weaker the shock 
wave (\ — m1), the smaller the o and the weaker the re- 
fection. Also o = | when A = 7/2. However, the 
flow behind the shock wave at AX = 7/2 (1.e., when the 
shock wave assumes the position of maximum flow de- 
flection) is subsonic, so that the concept of reflection is 
no longer a proper one. ; 

Finally, we should like to know the shape of the shock 
wave as a result of interaction with Mach waves. <Ac- 
cording to Eq. (45), 


¥+ 1 l 6p 


Say = - (45) 
2 sin 2a piyM? 
though admittedly an artificial 


A physical meaning 
Considering the 


one-—can be ascribed to the relation. 
pressure 6p/p,y.A/,° as the stress and the deviation of 
the shock angle from its normal position dw) as strain, 
then Eq. (45) can be interpreted as some sort of a 
with the modulus of elasticity 
Such a concept is useful in pre 


stress-strain relation 


2/(y + 
dicting the shape of the shock wave when the pressure 
immediately behind the shock wave is known and vice 
In many cases, the pressure immediately be- 


1)] sin 2a. 


verSa. 

hind the shock wave can be predicted quite simply; 

for, according to Eq. (21), the pressure immediately 

behind the shock wave is given by 

( 6p ) | 

piyM, at (2, x tan wi) \ M,;? — | 
6, tan w) | - G{x(1 + pb, tan wi) ]} 


}F[x(l — 


Using the boundary condition, Eq. (48), at the shock 


wave, we have 


cea) 
hyM,? at (x, x tan w1) 2 





(1 + o) tan wFl[x(1 — B, tan a) (56) 
so that 
‘ y+ 1(1+ o) tan p 
raat (2, = tam on) ™ 2 sin 2, 
F{(x(1 — B, tan w,)}] (57) 


Thus, the local 
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shock angle and the pressure immediately behind the 
shock wave at a point P on the shock can be related 
that is, the out- 


going wave F generated at the point Py: |[x(1 — 
This leads us to conclude 


directly to the incoming wave at P 


6, tan w,), 0] (see Fig. 4). 
at once that a compression Mach wave (generated at 
Py), when it interacts with the shock wave, will increase 
the shock angle whereas an expansion Mach wave will de- 
crease the shock angle. 

In fact, the above results are exact for a truly infini 
tesimal Mach wave, since the linearized theory of a 
steady two-dimensional flow involving only one Mach 
wave is actually exact in the limit (that is, the differ 
ence between the linearized solution and the exact 
solution can be made as small as one pleases by making 
the Mach wave weak enough). 

As another application of Eq. (57), let us consider 
the special case in which the reflected wave of the Mach 
wave generated at the tip of a small protuberance on a 
wedge profile hits the wedge again at a point down- 
stream of the trailing edge of the small protuberance, 
as shown in Fig. 6. In this case there are only out- 
going waves generated at the surface of the protuber- 
ance, and none of the reflected waves hits upon it. 
The small protuberance therefore lies completely in an 
irrotational field with outgoing waves only. It is well 
known from Ackeret’s theory that the outgoing disturb- 
ance F(x) is equal to the slope of the protuberance, so 
that the increase of the shock angle at a point P is 
directly proportional to the slope of the boundary at 
the point of intersection of the forward propagating 
Mach wave through P and the boundary, In other 
words, for this case the boundary will cast an exact 
‘“‘tmage’’ on the shock wave, except that the longitud1- 
nal scale is stretched by the factor sin (7 — w)/sin 
(u; — w,) and the transverse scale is changed, according 
to Eq. (57), by the factor 

y + 1 tan yp, sin (41 — @)) 
- (1 + o) = 
. sin 2w; sin py 


or 


¥y + | sin (441 —_ Ges 


2 sin 2w COS Ly 


(1 + o) 


In the general case the above conclusion does not hold, 
because reflection at the surface of the protuberance 
generally also enters the problem. However, since 
these reflections are generally weaker in amplitude, we 
may still say roughly that the boundary will cast a 
“slightly distorted image’’ on the shock wave. 

Since a compression Mach wave, on interacting with 
the shock wave, increases the shock angle, the entropy 
change for a particle is smaller when the particle passes 
through a compression. shock and a compression Mach 
wave than when the particle passes through a single 
stronger shock. This is also evident from Eq. (46). As 
a matter of fact, this is precisely a degeneration of the 
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principle upon which Oswatisch’s multishock super- 
sonic intake is based. On the other hand, by similar 
arguinent, one sees that the entropy change is greater 
when the flow undergoes first a compression shock 
followed by an expansion Mach wave than if it under- 
goes a single weaker shock. This statement is trivial 
once we have established that the interaction of ex 
pansion Mach waves and shock wave will reduce the 
shock angle and, hence, the strength of the shock wave. 


THE COMPLETE ANALYTICAL SOLUTION 


We have seen that the two equations governing the 
functions F and G are 


F(x) + G(x) = f’(x) (47) 


G(xl) — oF(x) 0 (48S) 


On eliminating G, we obtain a functional equation for 
F, 
[’(xT) (49) 


F(xT) oF (x) 


That is, we want to find the form of a function F(x) 
such that, when its argument x is replaced by xT’, the 
difference F(xT) — oF (x) is equal to a given function 
f’(x). The general solution of the above equation 
can be readily constructed (see Appendix). However, 
the particular solution corresponding to the present 
problem can be found without referring to the general 


solution. Such a solution is given by 


F(x) = f’(x) + of’ (x¥/T) + of'(x/T?) 4 


as may be verified by direct substitution. The deriva 
tion of this solution and the convergence of the series 


will be shown below. Rewrite Eqs. (47) and (48) as 


F(x) f"(x) oF (x/T) (58 


G(x) oF (x/T) (09) 


Let us suppose, for the time being, that /(.) is known 
This will, in fact, be the case if /’(x) 
xy, because this means that the protuber 


for x < %. 
O for x 
ance starts off at a finite distance from the tip instead of 
right on the tip so that /(v) = O for x < x. However, 
even if f’(x) # Ofor x < Xo, the following argument still 
applies.| Then, for x lying in the interval x) < x 
xol’, F(x) will be given by 
F(x "(x oF (x/T 
G(x oF (x/T 
is now known because of the argument 
r’, F(x) and 


where F(x /T 
2/T < 2%. Similarly, for «oI < x < Xo 


G(x) are given by 


: : ee: i 
F(x) Ag ® t ol ° t 0 fk ( E 
I I 
: De we 
G(x al ( ) o°F ( 7) 
I I 
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in accordance with Eqs. (58), (59), and so on, 
Xol"t!, F(x) and G(x 


In gen- 


3 . aT? ~ ° . 
eral, for xoI” < x are given by 


r Ir" 


G(x) 


ef: 
ao} E 
I 


However, since F(v/I'" +!) is actually not known, thes 


cannot be counted as the “‘solutions,’’ except in sucl 


cases as /’(x) 0 for x < xp. Since o |, it is clear 
that the latter terms in the above expansions contribute 
less and less if (a) the slope of the protuberance ['() 
is bounded and (b) the function F(x) is bounded fo 
x <x. Now, for a fixed x (instead of a fixed xp), the 
condition xl" < x < xl” 
trarily large n by taking x» small enough so that Eqs 


(60) and (61) still give the values of F and G at that 


t! can be satisfied for arbi 


particular v. Intuitively, as nm ~ ©, the contributior 
of the last term o’ 


We should therefore expect that F 


R(y/T" +!) in these equations should 
be negligible. 
and G(x) can be given by the infinite series 


, ; v » , Xv 
F(x) P(x) + of ( ) { ao-} ( ‘) T 
I I 


Gil Y 


To verify this, we have to establish first of all that 
these series are meaningful that is, that they are col 
vergent and therefore have a limit denoted by /; a 
G,; here. Next, we have to show that the differences 
F Fi and G 


assigned small number by making 7 large enough. 


G; can be made smaller than any pre 


If /’(x) is bounded for all x > 0, the series in Eqs 
(62) and (63) clearly converges both uniformly and ab 
solutely for o |, which is our case. Consequently, 
both series have a limit that has been denoted by F 

and Gi(x That F(x 
proaches infinity follows from the regularity requir 
iS bounded 


— F(x), G(x) > G(x) as n af 


ment at the tip of the wedge—namely, F(x 


in the neighborhood of x Q. 


= ‘ . - 9) th 
Since x) and » do not occur in Eqs. (62) and (65), ti 


last solution is applicable to all x's. 


known, we can calculate all the 
Thus, the 


With F(x) and G(x 
quantities that are of physical interest. 


pressure field is given by Eq. (21 
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Sh Y M,? rr i v =?) {tv hy 
a on ris Biv) tall ( 2 + / ( : ) 
pn Vue—it 
Spee: a) (: f nw) qf A) (x4 4 ) 
( / + + ¢ f + / Bra ie (64) 
| ( r : r a i, rs ‘ar ; p * 


In particular, on the wedge surface, the pressure distribution is given by 


6p Ei (*) r(*) (=) ; ; 
: I (et 4 2oa| ‘ t 2o°/ = r 2o° : ree (69) 
piyM, VM,2-1L I I I 


\sa check, consider a point x x, so near the tip that /(w) = /’(O) for alla < x; Then, Eq. (65) becomes 
6p/piyM? f'(0)/V M,? 1] (1 + 20 + 20? + 2o 
hat 1s, 
6p/ pry M, f’'(0)/V M, l Ll + a)/( 0 


When the value of o is substituted into the above equation, we obtain 
6p/piyM, /’(O) tan Xv 


which agrees with Eq. (43). 
Let us next estimate the effect that a small protuberance on the wedge has on pressure distribution on the wedge 


downstream of the protuberayée. Suppose that the small protuberance terminates at a point x; then, 


for x1 * o< Sit "(x 0 
ora < 2 < xt", f"(x ris/t 0 
for 40? < x < mF, "(x "(x /T (*(x/f 0 
w 
for xI* < x < x,J*t, "(x f"(x/T 47 f(x /T? (0) 


Hence, for a point x > x1", the pressure at the wedge surface 1s at most of the order 


bp A ll ie(=.) p (=) y’ (5) | 
piyM, VM peti ; pete ( pets iran 


that is, at most of the order of o4*!. The direction of the stream line at any point is given by Eq. (20). 


ye ; : (" mr) , (; t *»)) | 7 (" ) ? (; t *)| = 
P(A BY) -- e17 / t oO rs. SI 
My ; 4 r J : L : 


lhe shock wave angle is given by Eq. (57 
(y+ 1) sinA sing [ . ( sin (yy ray ! (* sin (uy w ) v sin (uy “)) | 
buy =~ f' (x — + o aor bef t— — mas 
sin 2w; sin (A + mw) LU SIN uy COS wy I’ sin wy COS w I sin yy COS a iH 
(67) 


because 


sin by W@W) 2 sin A cos Mi 
v(1 6, tan a, oo , r 6 : 
S11] 4 COS w@ sin (A + wy 


rhe entropy distribution is a function of yonly. Using Eqs. (46) and (56), we find 
6S y(1 M,? sin? @,)* sin \ sin yy f? ( sin (py) Ww} ) . (° sin (yy; Ww) ) 
= : y : : + al : : ; 
4 | eet sin (A + wa) Lo " sin py SIN @, [' sin wy Sin w 
W7,° sin? a 


R !' 
sin- w, {| 1 4 
») 
, fy sin (py ar ; 
ot : - ee oe (HS) 
I S101 wy SIN @) 
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Flow over an airfoil at an angle of attack. 


Comparing Eqs. (22), (23), and (24) with Eq. (21), we 


see at once that the density, the temperature field, and 
the variation of 64/1, are given by 


Pi piyM-? c p 
6T a 6p 6S a 
Tr, = (y — 1)M;/? i 4 C, (70) 
bu = 6p 7 I (=) (71) 
My piyM? yMir\R 


SUPERSONIC FLOW Past A THIN AIRFOIL AT LARGE 


ANGLE OF ATTACK 


A practical application of the above analysis is the 
case of a supersonic stream past a thin airfoil at a large 
angle of attack, as shown in Fig. 7. 

When a thin airfoil is placed in a supersonic stream 
at a large angle of attack, the flow at the top surface 
will, in general, consist of an expansion fan at the lead- 
ing edge as shown nr Fig. 7. The flow then smoothly 
follows the airfoil surface until it feels the effect of the 
trailing edge at which the shock wave is produced. 


SCIENCES JUS, 


In this region the flow is isentropic and irrotational 
Above all, only one family of characteristics jg jy 
volved, so that the exact analytical solution of the floy 
can be easily constructed even without introducing 
linearization. Alternatively, a solution can be easily 
written by combining the Prandtl-Meyer flow and th 
Ackeret’s solution. For flow that 
airfoil, a shock wave will be produced at the tip of the 


passes under th 
wedge. This shock wave will interact with the Mach 
wave generated at the lower surface so that the flow 
When the shape 
of the airfoil does not deviate too much from a straight 


becomes rotational and anisentropic. 


line, our analysis provides the analytical solution of g 
region of flow between the leading-edge shock way 
and the first Mach 
analysis, the resultant drag on the airfoil and the differ 


trailing-edge wave. From th 
ential pressure force on the upper and lower side of the 
airfoil can readily be calculated. 

As we have seen, the pressure distribution on the wall 
As o are usually small, we shall 


Thus, 


is given by Eq. (65). 
retain only the first two terms of the formula. 


5 Hy | + 
* : i } zal’ (*) | 7 
piyM, VM? — 1 


The first term corresponds to the pressure distributior 
obtained by Ackeret’s theory. The second term repr 
sents the effect of reflection from the shock wave on the 
pressure distribution on the wedge. In Fig. 8, thes 
effects are drawn separately in (a) and (b) for a para 
bolic are airfoil. It is clear that, for such an airfoil 
(or a similarly shaped one), there will always be an in 
crease of wave drag and a forward shift of center oi 
pressure. 

When the angle of attack is not very large so that the 
shock wave is weak or moderate, the pressure distribu 
has been calculated by 


tion on a “supersonic airfoil’ 


Busemann and later by Donov.'® 


APPENDIX THE GENERAL SOLUTION FOR THE 


FUNCTIONAL EQUATION 


We have seen that the equation governing the fun¢ 


tion Fis 
F(xT) — oF(x) = f'(x) (A-1 


This is a linear nonhomogeneous functional equation 


for F. Let us suppose that F = Fi, F = Fy are two 
solutions of Eq. (A-1)—.e., 
F\(xV) — oFi(x) = f'(x) 
Fi(xT) — o F(x) = f'(x) 
Then 
Fi(aT) — Fo(xl)] — of Fi(x) — Fi(x)] = 0 


That is, the two solutions F;, F. can differ at most by 4 


solution of the homogeneous functional equation 


F(xT) — oF(x) = 0 (A-2 
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INTERACTION OF 


Consequently, if F,(x) is the general solution of Eq. 
(4-2) and F,(x) 1s any particular solution satisfying 
Eq. (A-1), the general solution of Eq. (A-1) is given by 


F(x) = F(x) + F(x) (A-3) 


Before finding the general solution /,(x) of Eq. (A-2), 
let us first determine some particular solution of it. As 


4 trial, let us assume 
F(x) yl 


then Eq. (A-2) becomes 


( rr ym ox” 0 
[his equality can be satisfied for a// values of x if we 


choose m such that 


m log a/log I 


~») 


Consequently, a particular solution of Eq.(A-2) 1s 


. log o/log I 
F(x) i ati (A-4) 
A more general solution can readily be constructed if 
we let 


’ log a log I - 
Cis) = A-5) 


F(x) 
for, if we substitute Eq. (A-5) in Eq. (A-2), we obtain 


C(xP) C(x) A-6 
Hence, any function C(x) defined in the semiopen in 
terval |x, xo!) is a solution, provided that outside the 
interval the function is defined by Eq. (A-6). 

We shall now show that F(x) as defined by Eq. (A-5) 
and (A-6) is, in fact, the general solution of Eq. (A-2 
F(x). It is clear from Eq. (A-2) that, once F(x) 1s 
specified at x», then its values at xoI’, xoI"*, wT’, . . . are 
known. Now, if F(x) is specified for all x in the half 
opened interval |x», xl"), then F(x) will be defined for 
ill values of x > O or x Q depending on whether 
%) > O or Xy < 
allowed 


0; no more degrees of arbitrariness are 
Hence, if we restrict ourselves to the half 
line x 0, the general solution of Eq. (A-2) must con 
tain an arbitrary element such that the function F(x 
may take on any arbitrarily assigned initial values in the 
half-opened interval [x», xl"). Now, C(x) is-an arbi 
trary function which will serve this purpose; therefore, 
Eq. (A-5) and (A-6) which actually represents the gen 
eral solution of Eq. (A-2). Consequently, the most 
general solution of the functional equation, Eq. (A-1 


1s 
F(x) = C(x)x'® 78 4 Fix A-7 


where F(x) is a particular solution of Eq. (A-1 
For our problem, we have found a particular solution 


, 


F(x) = f(x) + of (x/T) — of (x/T?) +... A-S 


which remains finite at x = 0. Since o < 1, > 1, 
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Pressure distribution according to Ackeret’s theory. (b) Addi 


tional pressure forces due to reflection from the shock wave 
(c) Pressure distribution contributing to the force F,, perpendicu 
lar to the chord. (d) Pressure distribution contributing to the 
force F,; along the chord Note that drag F 


sin a + F, cos @ 


the exponent log ao log T° is negative. Hence, for 
F(x) to be finite at x 0, we must have C(x (0. 
Consequently, the only solution that remains finite 
at x = O (the fifth boundary condition mentioned in 
is F(x) = Fy(x is giver by Eq. 


the text , where Fy(x 


(A-S). 
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Downwash Behind a Two-Dimensional 
Wing Oscillating in Plunging Motion 


E. LAPIN,* R. CROOKSHANKS,? anv H. F. HUNTER? 


Douglas Atreraft Company, Ine. 


SUMMARY 


4 solution in closed form is found for the downwash behind a 
two-dimensional wing oscillating sinusoidally in plunging motion. 
The polar representation of the downwash at a given distance 
behind the wing is a divergent spiral. At large distances behind 
the wing the downwash variation with distance is sinusoidal. 


NOTATION 


wing plunging velocity 


h = amplitude of wing plunging velocity oscillation 
w = amplitude of downwash oscillation 
U = free stream velocity 


= airfoil chord 
Ky, Ki = modified Bessel functions 
= circular frequency of wing plunging oscillation 


i = 
t = distance in the wake 

Y = distance along the airfoil 
7 vorticity intensity 

H, = Schwartz function 

Hy, H; = Hankel functions 

Ci, Si. = cosine and sine integrals 


DISCUSSION 
| SHE BEHAVIOR OF THE DOWNWASH behind oscillating 
wings is of some importance in the study of air- 
craft dynamic stability. The downwash behind a two 
dimensional oscillating wing represents a limiting case 
in which only the influence of vorticity shed parallel 
to the airfoil trailing edge exists. 
von Karman and Sears! have shown that, for an air- 
foil oscillating with plunging velocity, 
h = Agl exp (ivt) (1 
there is a vorticity distribution in the wake 
—2rAoU exp (tvt — ivt/U) 


TS = Ewer poe 2 
Ko(iv l > + A \(av l ) 


and a vorticity distribution over the airfoil chord in- 


duced by the wake vorticity 


—2AoU exp (ivt) I — x 
= —— a coe . | 
| K(iv/U) + Ky(iv/l 1+ x 


“= exp (—ivt/U) le a 
= dé (3) 
JI ok inn 


There is also the quasi-steady vorticity distribution 
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Fic. 2. Polar diagram for downwash at a distance of five half- 
chords behind wing mid-chord. 


over the airfoil chord, 


V(x) = 2A UV (1 — x)/(1 + x) exp (zt) (4) 


By the Biot-Savart rule, the downwash at a point &, 


measured in wing half-chords, behind the wing mid- 
chord and in the plane of the wing, is 
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Str 
Str 


‘ 1 f° vdx Ly (' volx) dx “l vi(x) dx "= y() dé 
Wi &1) = = = 9 = is % 7 e { 
am. I aod: Qn |. 1 2 wt Zr =X Ji a = 
Substituting the expressions for yo(x), yi(x), and y(é), there results 
mw) ie (1 — x) dx | “ I can x {exp (—ivt/U) fé+ 1 
= SS 7 : a x 
h ig — x) V1 — x? Ao + Ai LJ-1 & — x Ee et '— x eo 
"© exp (—ivé/L’) dé 
didx +7 / | = ( 
J! le 


The first integral on the right, the airfoil quasi-steady vorticity distribution contribution to the downwash, js 


"" (1 = d 
/ x Xx =mril — V(& — 1)/(& 4+ 1 7 
1(& — x) Vl — x? 


The second integral, the contribution to the downwash produced by the wake induced vorticity distribution over 
the wing, may, by interchange of the order of integration, be put in the form 


o | I — x *? exp(—ivt/ LU’) le + | “= exp(—ivé/ LU’) le + 1&— 1 
= a & didxy =m - z ‘ee 1} dé (8 
FAG yt T+ ew [= 77 J\ &—- § e- tht ! 


It is seen that the second term on the right just cancels the downwash produced by the vorticity distribution in th 


wake. Thus, 


w(é;) i | | “> exp (—ivt/U) fe + 1 
— ‘fom = i+ : I ° dé () 
h & + 1 Ko + Ay Ji & —¢ &é—] 


Denote the integral on the right of Eq. (9) by F(z, £1), where s = v/l’, so that 


ae "= exp (—isé) le +1. "exp (—isk) dé | ; , exp (—ist) dé 
F(z, &) = - : - dé = — + (1 + &) 10 
J! a ¢ p= { ws “et — | a &=-ojVe | 


But 
exp (— vst) dé 
/ ——— Ko(iz 
JI Ve -—1 
and 
; — F exp (—Jst) dé 
F(s, &1) —Ay(iz) + (1 + &)J, Is, &) = J | ———— 
rf hove I 


The integral / may also be written 
, _ "= exp |—is(& — &)] dé 
I(s, &:) = exp (—/2é,) (1? 
di m=- 27 =— I 
from which 
dl . f° exp (—Isk) dé re - , 
= —1i,J+1 / = —16,/ + 1Ao(tz) = —7&/ 4 Ha’ (s 13 
dz J Ve -— 1 2 
A solution to Eq. (13), with exp (7&3) used as an integrating factor, is 


, 5 exp (—7s&,) dé he a 1 3 
I(z, &) = = /(&) exp (—7813) + exp (—7813) exp (& 0) L107 (0) dé 14 
i (f, — £) Y £2 — |] 2 0 


sl Ss s 


The constant of integration, from Eq. (11), is 


3 dé In (& + Ve? — 1) ; 
HO, &) = f&) = | = 1 o1 13 
D tee) Vee — I Vi? — 1 
and 
In (& + Vé&? — 1) xr {? ; 
I(z, &:) = exp (—?&3) ° = Ts | Hy"'(¢) exp (ie) dé (16 
Vit — 1 2 So 
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gy using the transformation £¢ = w, the integral appearing within Eq. (16) may be identified with the functions 
‘vestigated and tabulated by Schwartz.” 
T "* a). ie T a ae a 
Hy’ (€) exp (t&) dg ~ 7, a ', of) (17 
2 Jo 26 
The downwash to angle of attack ratio is then 
; ; In (£; + V&,"? — 1) T * 
alt, _ | —Ko(iz) + (1 + &) exp (— 7&2 - = +— &, ("68 
— om |” | Ve" - 1 2E1 (IS 
h G1 1 + = = : “ | 
| Ko(tz) + Ay(iz) 
But Ko(?z ~(im 2)H,) (zs) and K,(iz) = —(m/2)H,(s), so that 
w(E ws i te le, i” le x 
j Y: a \:, fe. 
a teat VWeP—1). t co 
(Mo (s) + (1 + seo : = + - H.(&-, &) 
T Vie - 1 £1 (19 


H,°)(z) + tH (z) 


In Fig. | there is presented the variation of downwash with distance aft of the wing mid-chord for several values 
of the reduced frequency parameter (wc/2l’) = s. At zero reduced frequency, Eq. (19) reduces to 


@(f)/h\,-0 = 1 -— V(ii — 1) ge 2) 


which is the downwash due to the quasi-steady vorticity distribution alone. A polar representation of the down 
wash is shown in Fig. 2 for the particular location £& = 5.0. The polar is a divergent spiral that may be demon 


strated from Eq. (19) by using the asymptotic values for H7/,°’, Hy)’, and H,;’. Schwartz gives 


9 2 | 
H,(&~', ©) = 3 E + In (& — V a 1)| 
Vi? —1 7 


which, together with 


Hy? @ (Cl + 1)/V 2s] exp (—iz), My” @& [(—1 + 1)/V 22] exp (is) 


yields 
oe me | 
w(&, | f&— | TS an : 
Pe ~~ s ™ e i[2(éi—1 w/4)] (20 
. i 2 We +1 2 
Eq. (20) indicates that the downwash becomes infinite as the square root of the reduced frequency. This is illus 
trated by Fig. 3, in which the downwash reduced by the factor zs” “is plotted. At extremely large reduced fre 


quencies, the diagram becomes a circle with radius (7/2) *. The physical reason for the occurrence of infinite 


downwash is that the vorticity deposited in the wake is proportional to dl’ /dt and dI’/dt is proportional to the fre- 
The infinite vortices in the wake produce infinite downwash. 


quency |see reference 1, Eqs. (21) and (23) |. 
1 


For calculating @/h at frequencies outside the range of the tabulated function H, (&,~", & 2), the following 


asymptotic expression of this function, given by Schwartz, is useful: 


E (kt } p—l ! + 
) (1 — 1)& exp [7(& — 1)z] —1 
H.?(¢,—', &2) = H.?(&7!, @ ) + At 4 R, | 
(£; — 1) V x2 »=0 i2(€; — 1) 
A, = [(2v)!/v!] D5 [(2m)!/(m!)?] [— (& — 1)/8]” 
m=(0 
R, = O(1/2°) 


The behavior of the downwash at large distances from the wing is of some interest. Returning to Eq. (6), for 


§: > x the second integral is approximately 
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—] a 1 — “? exp (—iszt) — /E | 
(RK / ee — dé dx 
&(Ko + Aj) l b= x =x =i 
iil . - le ra | >] r 
— : exp (—7z&) - dy 
(Ko + Ay) S1 gé- 1] 1&€ — x 1+ x 
—1 its ~é+ 1 jé — 1 
= —_—-; : exp (—72é) t— dg 
£1(Ko + Ky) l é — | é + | 
From reference 1, this is seen to be 
= : : exp (— 72 
a = = Ky Tr A, — : 9 
£(Ko + Kk) 1Z 


The third integral in Eq. (6), taking the Cauchy principal value, is 


—T “= exp (—izt) dé —rexp (—izh)[ £9 oles: T 
: : * = : : Ci2(&; — 1) +7281 32(& — 1) +2 29 
Ko + Ky 1 & — = Ko + kK, ys 


Thus, the downwash at large distances from the wing is approximately 


Ww 


¢ e 1zZ 
os Sm = cee 
h | | & | 12(Ko + =| 


and for &, very large 


ctr 
| 


utr 
+ 


_ . at 
v ime **' 
~ (24) 


h Ko(iz) + Ky(iz) 


a) 


The polar diagram, as a function of the distance &, 
also tends to a circle. This is illustrated by the dotted 
curves in Fig. 4. The physical interpretation of the 
constant frequency trace is that at large distances be- 
hind the wing the downwash alternates sinusoidally in 
space, as well as with time. 

Lingard’ has also calculated the downwash behind 
an oscillating two-dimensional wing under the approxi- 
mation that the bound vorticity on the wing is concen- 
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illustrating asymptotic behavior as wo/2U — 


e 1261 C; x a rT 
a> : Be) — 1) + 1 Sh - ie ee. 23 
Ko + A, e F 4 


trated at the mid-chord. His results and these are in 
good agreement. Nocomputing advantage is gained by 
the approximation, since his final formulas are not 
simpler than those given here. 

Calculations are in progress on the downwash behind 
oscillating finite aspect ratio wings using an extension 
of Reissner’s' procedure for the calculation of the load 
distribution on such wings. 
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_ Analysis of Elastic Structures by Matrix 


> 


Transformation with Special Regard to 
Semimonocoque Structures 


BORJE LANGEFORS* 
Saab Aircraft Company 


SUMMARY 


The system of equations for a redundant elastic structure is es 
tablished by The is divided into 
separate steps, the first being the establishing of the ‘ransforma 
This first step is the only 


a routine procedure work 
sion matrix by purely statical analysis. 
one where knowledge of structure analysis is demanded, all fur 
ther steps being simple matrix multiplications that can be per 
formed by automatic means such as punched-card machines 
The procedure is a modification of G. Kron’s method for electric 
networks 

The structure can be analyzed by parts. Subsequent changes 
in member stiffnesses do not affect the transformation matrix, 
which can therefore be established once the geometry of the struc 
ture is specified. The rest of the work is a simple routine 
matter 

The system of equations obtained is reduced by matrix trans 
formations of a similar kind. The procedure for obtaining the 
system of equations is illustrated by a simple example 


INTRODUCTION 


: be RE IS A CERTAIN RESEMBLANCE between the 
analysis of an elastic structure and that of an elec 
trical network. In both cases simple members are 
coupled together to constitute more or less complex 
systems. The problem of analysis is in both cases that 
of finding the physical state of each element, this state 
being a consequence of the introduction of certain dis 
turbances into some parts of the structure. The key 
to the solution is in both cases the manner in which the 
constraints, caused by the coupling together of the 
influences the distribution of the 
The most general types of 


different members, 
disturbances introduced. 
disturbances are for both cases of a dynamical character 
and such dynamical problems with constraints are 
solved in a general way by the method of Lagrange, 
using his generalized coordinates. This method was 
introduced by Lagrange for mechanical systems and 
was used by Maxwell in a modified form for electrical 
systems also. It may be interpreted as a coordinate 
transformation (see reference 1, pages 139-140). In 
the static cases encountered in elastic structure anal- 
ysis and in stationary network theory, the method of 
Lagrange simplifies to the methods of minimum poten 
tial energy or minimum work. This is seen immedi- 
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ately from the form of the Lagrange equation (see also 
references | and 2). 

More recently, G. Kron 
method for the analysis of electric networks. 
method could be interpreted as the introduction of the 
Lagrange generalized coordinates by means of a com- 
According to 


another 
This 


has introduced 


mon transformation of coordinates. 
Kron, this transformation is defined by the conditions 
of constraint put into a matrix relation between the 
primary coordinates and the generalized coordinates. 
Because of the close relation between the Kron and the 
Lagrange methods, it is to be expected that the Kron 
method will be applicable to elastic-structure analysis 
and to dynamical problems. 

FOR ELECTRIC 


TRANSFORMATION METHOD OF KRON 


NETWORKS 


The method of Kron for electric networks consists, in 
short, of writing down the equations of performance for 
each element, these being expressions of Ohm's law. 
When these equations are written in matrix form to con 
stitute the “‘primary equation,” the currents and the 
voltages related to each element will occur in one cur 
rent- and one voltage-column matrix, and the resist 
ances or conductances are placed as elements of the 
principal diagonal of a square matrix, the other ele 
ments of which are zero. Then the relations between the 
currents or voltages in each element and some chosen 
real or hypothetical currents or’ voltages, serving as 
generalized coordinates, are written down in matrix 
In this way the matrix of the ‘‘coordinate trans- 
Then the rest of the work to 


form. 
formation” is defined. 
find the equations for the network consists of a routine 
procedure using the rules for multiplication of matrices. 
Kron has presented his method in references 3 and 4. 
The author regards reference 7 as the best introduction. 
(See also reference 8.) 

THE KRON METHOD TO THE ELASTIC 
STRUCTURES 


TRANSLATION OF 


In order to follow Kron’s procedure in the elastic case, 
we first must write the ‘‘primary equation’ expressing 
the physical laws for each elastic element. It will be 
found that for some elements the physical laws will be of 
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exactly the same form as in the electrical case, whereas 
for other elastic parts this will not hold. In the latter 
cases, however, the same types of equations will still 
be found if these are written in matrix form rather than 
in scalar form, and simple letter symbols are used for 
each such elementary matrix. These elementary ma- 
trix equations will constitute some generalized form of 
Hooke’s law: 

vP =u (1) 


v, P, and u being scalars or matrices, indicating flexibil- 
ity, force, and displacement, respectively. 

The next step in the analysis will be to establish, in 
matrix form, the relation between the forces P; acting 
on each element and the redundancies X,, the latter 
being the generalized coordinates. In that way a re- 
lation of the form 


P; = QiaxX) + Aik. + 


is obtained which, in matrix notation, becomes (com- 
pare reference 7, page 13, or reference 11, pages 6 
and 8) 


P = CX (2) 


P and X being, respectively, the column matrix of the 
forces in all the members of the structure and the col- 
umn matrix of the redundants. (For the simple ma- 
trix operations used in this paper see references 4, 7, 
8, or 11.) Cis the transformation matrix of coordinate 
transformation. 

Concerning the transformation (2), it is to be ob- 
served that C is a rectangular matrix and, hence, is 
singular—i.e., it has no inverse. Therefore, there is 
no one-to-one correspondence between the 
nates’ P; and the ‘‘coordinates’’ XY; as would be the 
case were Eq. (2) a real coordinate transformation. 
However, it is possible to add to the Y; so many new 
components X, that the number of X; + X, would be 
the same as the number of P; and in such a manner 


“‘coordi- 


that the quadratic matrix so obtained would be non- 
singular. It is the am of introducing the generalized 
coordinates of Lagrange to get a reference frame where 
only so many of the coordinates are different from 
zero as there are degrees of freedom—4.e., for our case 
as many as there are redundants. 
show that the use of the transformation matrix C will in 
all operations used here give the same result as would 
the use of the nonsingular matrix, together with the 


It is then easy to 


conditions XY; = 0, which are just the conditions of con- 
straint of the structure. The use of C will have the 
advantages that the calculations become shorter and 
that it will not be necessary to define the variables X 
exceeding the redundants_X ;. 

The ‘‘primary equation” will be, in matrix form, 


VP = 4 


where | indicates the matrix whose diagonal elements 


are the scalars or matrices v of the simple members oj 
the structure and where wu is the column of aj the 
displacements or displacement matrices of the mem, 
bers. 


If the assumption is made that the Kron theory ap- 
plies to the elastic problem, as will be proved below, the 
system of equations for the structure will be 


C'VCX = C'u 


or 


if we put C’VC = V,andif C’u = u,. C’ is the trans 
pose of C (see reference 11, page 3). 

The validity of the Kron method for the elastic cas, 
remains to be proved. As Kron’s presentation of his 
methods does not seem to constitute rigorous proof, the 
demonstration will be given here in a somewhat differ. 
ent manner. This will at the same time give the ex 
pression for the u-column, which in the elastic case js 
seldom given directly, as in the electric networks, but 
must be calculated from the given exterior loads. Fur 
ther, the equations of state of the simple elements are to 
be given. This will be done here only for the elements 
that are encountered in the design of semimonocogue 
bars fastened to shear 


structures—i.e., pinjointed 


carrying plates—and rectangular and trapezoidal plates 
For other members see reference 10. 

It may be mentioned for completeness that Kron has 
established electrical networks, equivalent to elastic 
beams, and has pointed out that his transformation pro 
cedure for electric networks may be used for the anal- 
ysis of such elastic structures as can be built up by 
such beams.® The present author has later on estab- 
lished electric networks analogous to shear-carrying 
plates,’ and, consequently, we have at our disposal all 
that is needed to use the common Kron procedure for 
electric networks to analyze the networks that are analo 
gous to semimonocoque structures. These analog 
networks are, however, established by the aid of stiff- 
ness matrices, instead of flexibility matrices, which will 
lead to a much higher number of unknowns. Further, 
in this case, the transformation from the primary matrix 
to the final one is so simple that no matiix transforma 
Whereas the stiff- 


ness matrices seem to be the best means to establish 


tions are needed for its execution. 


network computers for elastic problems, the direct 
calculation, by digital means, is better performed by the 
use of the flexibility characteristics of the structure as 1s 
the most common way in redundant structure analysis. 
This cannot be performed in this case by simply invert 
ing the stiffness matrices to give flexibility matrices, 
because the stiffness matrices mentioned are singular 
and so cannot be inverted. Hence, it is necessary 
to follow different paths, as has been done in this 


paper. 
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ANALYStis OF 


ProoF OF THE TRANSFORMATION METHOD APPLIED TO 


ELASTIC STRUCTURES 

For a simple elastic element, where Eq. (1) is scalar, 
the load P will, when applied to the element, do the 
elastic work (1/2)Pu, which, according to Eq. (1), can 
be expressed in the form (1/2)P?v or (1/2)PvP. This 
work will be equal to the elastic energy stored in the 
element. As will be shown later, in connection with 
the establishment of the equations for the elastic ele- 
ments, the energy UL’ in the parts that are guided by ma- 
trix instead of scalar equations is also given by an equa- 
tion of similar form, so that for all elastic elements we 
have 


2U = PvP 


where P and v are matrices or scalars, depending on the 
type of the element. P’ is the transpose of the matrix 
P. 

For the whole elastic structure, the elastic energy 


will, of course, be equal to 
OU; = (1/2) PoP 


where the index 7 denotes the successive numbers of the 

elastic parts of the structure. It is seen that the total 

energy is a quadratic form and can be expressed in ma- 

trix form as 

- a die 
2 
P, 


Do P, 


a 


P, and 2, still being scalars or matrices depending on 


2U = [P,'P.’ ... (4a) 














= — 


tvpe of element), or, in a condensed form, 


2U = P’VP (4) 


see reference 11, page 28). 

When the structure is of redundancy , it is neces 
sary to introduce k unknown forces, the redundants X, 
in order to calculate the compound column P, the ele 
ments of which are forces or force columns P;. Then 
it is possible to express the forces acting on the elements 
of the structure as linear functions of the exterior loads 
Q,; and the redundants X,, this relation therefore being 


in the form (to be compared with Kron’s Eq. (2) 


? ) ‘’D * 
I CX + DQ ce O | 


CD] and P,.’ = aa 


or 
P = C,P. where C; 


Now, obviously, the energy L’ can be expressed in terms 


of the loads and the redundancies-~i.e., by P,and, as 
it must still be a quadratic form, in matrix notation it 


will be 


2U = P.'WP. (4c) 
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where IV remains to be determined. From the equality 
of the expressions for the energy, we get 
P.’WP, = P’VP = P.'C,' VCP. 
[from Eq. (5)] and so 
W= C/VC (6 
The equations of state of the structure are obtained by 
the condition that the energy is to be a minimum when 


the values of X; are such as to maintain equilibrium. 


This condition will give the equations 


0U/oX; = 0 


a Cre CVD Tt x 
)\ | , ‘) d) 
(1/2) [X’Q paeaes | (4d 


we thus get the set of equations 


With 


U = 


OU/OX = C’VCX + (1/2) (C’'VDQ + Q'D'I'C) = 0 


(compare reference 8, page 48, where the notation 
{0/OX;{A is used instead of 0A/OX of the present 
paper). 

The two expressions in parentheses on the right-hand 
side of the last equation are seen to represent the same 
bilinear form (reference 11, page 28) and are therefore 


equal, and so we get 


OU/aX = V,X + C’VDQ = 0 (7 


where 

C'vVc (S) 
Now DQ is the column of the forces on the different 
members induced by the exterior load column Q, and 
thus DQ represents the displacements of the ele 
i.e., the column — wu (with the minus sign chosen 
Hence, the equation can be written 


ments 
for convenience). 


as 
VX = u, (9 


if we use the symbol u, = C’u. uu, is thus the column 
matrix of the displacements related to the redundants 
X, induced by the given exterior loads Q;. Eq. (9) is 
seen to be identical with the Kron Eq. (3), which is 
thus proved to be in accordance with the theorem of 
minimum energy (and, hence, with the theory of La 
grange in its simple special form for static systems). 

In the case of thermal expansions of the members, 
the corresponding “‘primary”’ displacements u, of the 
members are obtained at once from the temperature 
changes and thermal expansion factors of the material. 
Then the column matrix u, of Eq. (9) is again obtained 


as u, = C’u. 


PRACTICAL PROCEDURE 


The procedure for routine calculations will be the 
following (with some slight deviations from Kron elec 
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trical network procedure): The different steps (A), (B), 
and (C) to be discussed below are illustrated in detail 
by a simple example at the end of the article. 


(A) Establish the Relation 
P= CX + DQ = CP. 


and thus C; = [C, D]. 

This can be performed in symbols once the type of 
the structure is determined and numerically as soon as 
the form of the structure is defined, regardless of size. 
The most practical way seems to be to introduce P, (i.e., 
X, and the Q,) into the structure one at a time, all 
other X,’s and Q;,’s being kept equal to zero. The forces 
on the different elements thus obtained, for example, for 
X, = 1 are the corresponding elements of the 7th col- 
umn of the matrix C}. 

This is the only part of the work where insight into 
the theory of structures is needed, and the rest of the 
work can be done by inexperienced labor. Once a 
sufficient stock of standard matrices C and D has been 
procured, all calculations can be done by inexperienced 
labor and, to a great part, by punched-card machines. 


(B) Establish the Primary Flexibility Matrix V (Rather 

Than the Primary Equation of Kron) 

This is accomplished by merely writing the flexibility 
factors or matrices as diagonal elements of a matrix. 
Also, V’ can be written in symbols as soon as the type 
of the structure is laid down, but the numerical deter- 
mination of |’ demands all the stiffnesses of the parts to 
be determined. 
thicknesses of plates will change the V' matrix but not 
Ci. 


Thus, changes in the areas of bars or 


(C) Perform the Matrix Multiplications 


C'V=B 
and 

V, = BC 

F, —BD 


to obtain the final Eq. (7) 
V,X = F,Q 


(D) Solve the Final Equation by Some Standard 
Procedure To Give X as a Function of Q 
In matrix notation this can be symbolized by 


X = V,'F,Q 


(see reference 8, page 22). It may often be advantage- 


ous to perform a variable transformation 


» = A'V,A 


XA = AX ] 
40 = A’F.0S 


J 
F 


to obtain a new system of equations which is easier to 
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solve. This will be briefly discussed below (see alg, 


reference 10). 


(E) Obtain the Forces P, 


Obtain the forces P; in all members of the structure by 
inserting the solution XY into relation (5). 


THE DEFORMATIONS 


By establishing and solving Eq. (9), we obtain jj 
that is necessary to determine the forces acting on the 
different members of the structure. However, for 
aeroelastic purposes, for example, it is often demanded 
to calculate the deformations also. For this problem, 
the expression for the energy will again be useful. We 
have indeed, if « denotes the column of the displace. 
ments u; of the points of application of the different 


loads Q;, 

2U = u'Q = Qu (10) 
Now, if Z is the flexibility matrix of the given struc. 
ture—1.e., the matrix whose elements are the displace. 
ments of the points of application of exterior loads 
when these loads are equal to unity, 


u = ZO 
and thus 
2U = 0'ZO 12 


Now, in Eq. (4d) we have an expression for L’ in terms 


of X and Q. This expression can be written as 


i X 
: -s KS > | 


| 


with notations the definitions of which are immediately 


obtained by comparison with Eq. (4d). From Eq. (7), 


X = —V,—'-V,9Q 


whence, 
2U = QV eR! Vi" VV," VieqgQQ — QO’ Vig’ Ve" Veg? - 

0’ V9 Vz~' Vie + Q'V 00 
Comparison with Eq. (12) shows that for Z we get the 
expression 
V9’ Vi" Veg 13 
which, together with Eq. (11), gives the expression for 
the deformations. 

SOLUTION BY PARTS 


It will often be advantageous to subdivide the struc- 
ture into parts and to analyze them part by part, pro- 
vided the solution for the whole structure can then be 
obtained in a simple way. 


of great advantage if a part of the structure is to be 
changed also with regard to the topographic constitu- 
Thus, for example, it would be possible, when 


tion. 


Such a subdivision will be 
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ANALYSIS OF 


studying sweptback wings with different angles of 
sweepback, to analyze the outer wing once and for all 
and then to combine it with different interior parts of 
the wing, the different parts having different angles of 
sweepback and perhaps also different forms in other 
respects. It is the aim of this section to show that such 
a subdivision is often possible and how it may be prac- 
tically performed. 

Suppose that the structure shown in Fig. | is to be 


analyzed by the aid of subdivision. One of the sub- 


structures or parts is to be the one indicated in Fig. | - 


by I. 

“We suppose that the structure is of such design that 
, part P may be isolated from the other two parts I 
and II of the structure by k hypothetical cuts, leaving 
only three longitudinal stiffeners in each of sections 
To these hypothetical cuts are related 
If the section a-a is introduced 


a-a and b-b. 
k redundant forces p. 
close to a transverse stiffener (a-a) in such manner that 
this stiffener belongs to part I, it is seen that all the 
redundants Y, interior to part I influence part I only. 
Analogously, the redundants y, are seen to influence part 
Ilonly, whereas the unknowns /; do influence both parts 


[and II. From this we conclude that relation (5), 
P=CX + DQ, will become 
Py GC Dy, X dD, 
F, = Dd, p + py O (14) 
Pir Diy Cir y | D» 
Vy 
J i 


Vir 
D;, is the matrix expression for the distribution of the 
D,, and Dy, are analogously de- 
fined. The use of the letter D rather than C is be- 
cause the forces /, are for all three parts of the structure 


forces p into part I). 


to be regarded as unknown exterior loads. 
The transformations C’I’C and C’VD are now seen 
to lead to the following equation [analogous to Eq. 


7)| for the whole structure: 


V. Vz, Xx] —F,Q 
i Vy, Voli? | ={ —F0 (15 
L Von Wy ) F,Q 
where 
V.= CVG: Va = GV 
\* = C,’V,D, 
V, = Dip’ ViDip + Dopo’ VopDop + Duy’ ViPn, 
F (Di,p’ ViD, + Dop’ VirpDp + Dip’ VurD2 
[v. Ci Vii; Vp Cul nDu, 
F, — Cr’ Vind: 


From Eq. (15) we get, for Y and Y, 


Vx = —V,pp — FQ\ 


(16) 
V9 -Vypp — F,Q4% 


Both these equations are established when analyzing 


ELASTIC 


STRUCTURES $55 








parts I and II in the usual manner (described in the 


previous sections of this paper) tf only care has been 
taken to separate the loads p from the given exterior 
loads Q when establishing the D-matrices. 

It has thus been shown that, if the parts I and II 
have been analyzed previously, the equations of state 
(15) for the whole structure may be written down at 
once if |’,, and D,, are also established. This is always 
a simple matter. 

If the equations for parts I and II have not only been 
established before but also have been solved by inver 
sion of the matrices |’, and J’,, we can go still further. 
We then have 
'F.O\ 
‘FON 


Vp — V, 
‘Visi — V, 


x= —VJ, 7 
, (id 

y= —), 
from which we get, by inserting in the p-group of equa 


tions in Eq. (15), 


Eq. (18) is a system of k equations from which the & un 
knowns p can be calculated. 

We have thus found that, in order to analyze the 
whole structure when both the substructures have pre- 
viously been analyzed and the matrices of their equa- 
tion systems have been inverted, we only have to de 
termine the matrices J’,, and D,,, perform some matrix 
multiplications to establish Eq. (18), and solve this 
system of equations for the k unknowns /,;. 


MANIPULATION OF THE EQUATIONS 


The systems of equations occurring in the analysts ot 
often and, 
It is therefore usual to 


semimonocoque structures are verv big 
hence, troublesome to solve. 
reduce the number of unknowns by some simplifying 
assumptions. Two groups of methods for such variable 
reduction are common in mathematical engineering. 
One of these methods consists of regarding the system 
as built up of ‘‘stiff’’ subsystems so that the number of 
unknowns is reduced to the number of subsystems 
chosen. The other group of methods is related to the 


names of Rayleigh and Ritz and many of their followers 
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In the case of semimonocoque structures, the first 
method is usually applied by taking same longitudinal 
their number is reduced 


stiffeners together so that 


whereas the second group of methods is applied by 
choosing some force groups by which the solution de- 
sired can be built up approximately (of course, also 
exactly if the number of force groups is equal to the 
number of unknowns, all force groups being linearly 
Both methods suffer from the draw- 
back that the errors entered by the variable reduction 
The present author has previously 


independent.'” 


are not known. 


shown that both groups of methods can be regarded as 
By performing 


different aspects of one single method. 


the reduction by means of matrix transformations of 
the same type as those demonstrated above (applied to 
the final equations instead of introducing them already 


in the elastic analysis), one obtains the important ad- 


vantage that the errors can be determined, and thus the 
Further 


solution can be refined by successive steps. 
changes in the reduction procedures (for example, 
insertion of new approximating force distributions used 
in the Rayleigh-Ritz method) can easily be performed 
by simple matrix multiplications. It is also found 
that new types of reduction can be introduced by the 
transformation method (references 14 and 10). The 
reduction by the transformation method will be illus- 
trated here by a simple example only. Suppose stringers 
are cut by eight hypothetical cuts, and thus eight un- 
known forces occur. By taking together stringers | 
and 2, 3 and 4, and so on, we have 2; = X; = Xo, 2. = X3 
= x4, and so on (zg; being some sort of generalized co- 
ordinates), or in matrix notation 


- r - 














which can be written x = Az. 

If the system of equations is V,x = u,, we obtain, by 
transformation by A, 

(A’V,A)s = A’u, 

which is a system of four equations only. The solu- 
tion z of this system gives a step curve for x. How- 
ever, we prefer to use a smooth curve, drawn to replace 
the step curve, as an approximate solution x! of the 
By entering x! into the equations, 
we obtain u,' = V,x'. We then solve the equation 
system V7 = u,;1 — u, in the same manner as before to 


system V,x = Uz. 


obtain corrections to the x!-solution, etc. 

If we had preferred to use the Ritz type of reduction, 
but still use the transformation method, each approxi- 
mating group of forces would define one column in the 
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A-matrix, the rest of the procedure being the same as 
above. 

It is, however, now possible to introduce new types of 
approximation, still using the same transformation by 
only changing the A-matrix. One simple vet powerfyl 
way 1s to specily 


= X), 22 = XN3, Vo = (1/2) (21 - 2 


etc., from which we get 


be BE ' 














X7 

Ls} L 7 
As this is an expression in matrix form of rectilinear 
interpolations in steps, it is natural to try other known 
interpolation types. This may be a better procedure 
for some problems, but very often the simple rectilinear 
interpolation is found to be the most practical one. 

If the same A-matrix is used for all the consecutive 
approximating steps, it is practical to determine the 
inverse of the matrix A’I’,Z. Then the iteration pro- 
cedure will be easy to perform. 


EXAMPLE 


In order to illustrate the practical procedure, let us 
study the simple example of the structure shown in 
Fig. 2. 

Q is a given external load acting on the structure 
Let X, denote two equal and opposite forces—one 
applied on the lower end of bar 3 and the other applied 
Similarly, let .Y. denote 
two equal and opposite forces applied at the lower end 


on the upper end of bar 4. 
of bar 4. X, and X», are supposed to be compression 
loads and are therefore taken as negative in the follow- 
ing analysis. 

(a) To find the transformation matrix C,; = {[C D), 
we have to determine the forces that are induced in the 
elements of the structure —i.e., the forces atting on each 
end of the different bars or along one edge of each plate 

when X,, X2, and Q are acting in turn (with numerical 
value = 1), the others being zero. 

Fig. 3 is an illustration for the case when ., = 
QO = 0, X2 = 0. 
places in Fig. 3 represent the different elements of the 


The values written down at various 


i.e., the X,-column, of the transforma- 


the first element in the .Y,-column 


first column 

tion matrix C; 
i.e., the element Co, belonging to the first row—and 
the first column in C is the force acting at end point 
‘1’ on bar number “‘O1”’ (see Fig. 2). From Fig. 3 
i.e., Con, = 0. In the 
same way we obtain, from Fig. 3, Cow: = |, Con 

The first two 


we find this value to be zero 


—1, and Cyn = —2a 2b, for example. 
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ANALYSIS OF 


index numbers refer to number of member, the third to 
end point of member, and the fourth to the number of the 
matrix column.) By drawing figures corresponding to 
Fig. 3 for X» and Q, we thus obtain the complete 
matrix Ci. 

(b) Step (b) consists in establishing the primary 
fexibility matrix 1. We first need to determine the 
fexibility factors or matrices of the components of the 
structure—i.e., the scalars or matrices 7, 0%, ..., in Eq. 
(4), In the example of Fig. 2, all plates are rectangular; 
hence, all shear flows are constant and so the forces in 
the bars vary linearly. We therefore obtain, for. the 
energy in a single bar with constant area, 


= [ a dx = [P\(2P; + Ps) + 
2 Jo EA 12KA i 
P.(P; + 2P2)| 
where 
= length of the bar 
A = area of the bar 
E = Young’s modulus of elasticity 


the end forces acting on the bar 


P, and P, = 

This is a well-known relation; compare, fcr example, 
reference 9, page 531 [Eq. (9)|. We translate it now 
into matrix language to 


» > p 2r r P, 
(] 2) | F iP» |? bd || 


r = 1/6EA 


U = 


which can be condensed to 


U = (1/2)P'vP 
P’ = [P, P2| 


2r r 
v = 
Se 


P being the column matrix of forces acting on the bar 
and v being the flexibility matrix of the bar. 

It is shown in reference 10 that, for bars that are 
connected to the edges of trapezoidal plates, the force 
can still be regarded as varying linearly if the end 
loads are determined by moment equilibrium condi- 
tions for the plate. This procedure introduces two 
errors that are such as to very nearly balance each 
other. Hence, the )’-matrix shown above will be ap 
plicable for all bars loaded by shear forces along the 
length and with constant area. For bars of varying 
area and bars with bending and torsion elasticity, also 
see reference 10. 

For rectangular plates, the internal energy is given 
by the well-known expression 


/ = {j 2)q?(a-b Gt) 


where a and b are side lengths of the plate, / is the plate 
thickness, and g is the shear flow. 

If P is the force acting along an edge of the plate, 
Then 


we have g = Pb if b is the length of the edge. 


we can write the expression for L’ as 
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2)P(a/6GH)P = (1/2)P0P 


v = a/bGi 


=a 


v being the now scalar factor of flexibility for the plate. 
For trapezoidal plates we get v = (a/bGt) (1/¢), where 
¢ is a correction factor the numerical value of which 
depends on the trapezoidal form of the plate.” 

Now we are ready to write down the |’”-matrix, the 
primary matrix, for the structure of Fig. 2. We assume 
all areas of the bars to be equal and the same for all 
plate thicknesses. 

We have to insert into the matrix of Eq. (4) for 2 
the flexibility matrix of bar number -|-i.e., 


b 2 I 
~~ == 
OLA | = 
In the example chosen we have vy, = v2 = ... = %% 
and v; = vg =... = 4 = a/GEA. For the plates we 


get the flexibilities v3 = v4 = ... = vj, = a/bGt. Now, 
by entering the matrices 7, to vj. as diagonal subma- 
trices and 2; to v5 as diagonal elements of a matrix, 
we get the primary flexibility matrix |’ of the struc- 
ture. 

(c) The equations of state are obtained by _per- 
forming the matrix multiplications of step (c) above. 
As a result of this procedure, we finally get the equation 


system 
j a® KAa) jb as tAa\ 
ij 3 x 7 _— Xo 
Vt eT par fT Va 7 202 7 20 * 
yb a® ~~ EAa) 
af. ) 
4° 262 ° 20Gth * 
jb a® KAa) jo a’ KAa\ i 
4 oe 2G *' 2 7 302 7 aaah? > 
jd a® ( 
4 Gar 
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which is found to be in accordance with the result oh 
tained by a conventional procedure. 
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Creep Buckling of Columns 


CHARLES LIBOVE* 
Langley Aeronautical Laboratory, N.A.C.AA. 


SUMMARY 


The creep of a slightly crooked H-section column carrying a 
constant load is studied theoretically. The material of the col- 
umn is characterized by a strain-time relationship, under con- 
stant uniaxial stress, of the form e = (o/E) + Ac®??t’, where e€ is 
the total strain, o is the constant stress, / is the time, and E, A, B, 
This form was selected because it 
75S-T6 aluminum alloy at 600°F. 
However, the anal- 


and K are material constants. 
applies to at least two alloys 
and a low-alloy steel at 800° and 1,100°F. 
ysis is intended for any material having creep properties of the 
same form and for which the material constants are known. A 
strain-time relationship under variable uniaxial stress, necessary 
for the column analysis, is formulated from the constant-stress 
properties with the aid of Shanley’s engineering hypotheses of 
creep 

The analysis leads to the conclusion that the lateral deflection 


approaches infinity—that is, the column collapses—in finite time. 
Results are given showing the maximum length of time the col- 
umn can support a given load before it collapses and the growth 


of stresses, strains, and deflections prior to collapse 


INTRODUCTION 


= HIGH TEMPERATURES can develop in the 
skin of aircraft traveling at supersonic speeds, 
material creep may have to be considered in the struc- 
tural design of such aircraft, especially those intended 
for sustained or repeated flights. 

The creep problem in aircraft has two aspects not 
generally encountered in other high-temperature equip- 
ment. First, a large part of the aircraft structure con- 
sists of compression elements, such as columns and 
plates, which are relatively slender or thin and which 
unavoidably have slight initial crookedness or eccen- 
tricity of loading. Creep in such elements will tend to 
produce continuously growing lateral deflections under 
constant load and may eventually lead to collapse or, 
as the collapse of a column is generally denoted, buck- 
ling. Second, the service lifetime of the aircraft may 
be short measured in hours or minutes instead of years. 
Hence, most of the significant creep may be of the pri- 
mary type, which is nonlinear for constant stress, rather 
than the simpler linear secondary type. 

In the present paper a problem is considered which 
incorporates both new aspects—that is, the continuous 
lateral deflection of compression elements and a non- 
linear type of creep law for the material. The com- 
pression element considered is an H-section pin-ended 
column with a slight initial sinusoidal crookedness in 


the plane of the web. A load below the Euler load 


Presented at the Structures Session, Twentieth Annual Meet 
ing, I.A.S., New York, January 28-31, February 1, 1952 
* Aeronautical Research Scientist 


is assumed to be rapidly applied and then kept con- 
stant, and the subsequent deformation due to creep is 
studied. (The rate of load application is assumed 
to be sufficiently high so that the creep occurring dur- 
ing the loading period is negligible, but not so high that 
dynamic effects have to be considered.) The material 
of the column is characterized by a nonlinear strain- 
time relationship, under a rapidly applied constant uni- 


axial stress, of the form 
e = (c/E) + Ae*’t* (1) 


where e is the strain, o 1s the stress, / is the time sub- 
sequent to load application, and /, A, B, and A are 
This form was selected because it 
75S-T6 


material constants. 
has been found to apply to at least two alloys 
aluminum alloy at 600°F.!' and a low-alloy steel at SOO 
and 1,100°F.? 

Two major assumptions are made in the analysis. 
The first is that a creep law for varying uniaxial stress 
may be formulated with sufficient accuracy from the 
constant-stress creep law by use of Shanley’s engineer- 
ing hypotheses of creep.* The second is the simplifying 
assumption that the deflected shape remains at all times 
sinusoidal. In consequence of the second assumption, 
the requirements of equilibrium, compatibility, and 
material behavior need be satisfied only at the mid- 
height of the column, and the solution obtained is there- 
fore a one-point collocation solution of the actual prob- 


lem. 
SYMBOLS 
A,B = material constants used in describing constant- 
stress creep curve 
b = distance between centers of flanges of column 
E = Young's modulus 
K = material constant used in describing constant-stress 
creep curve 
5 = length of column 
a = cross-sectional area of column 
t = time; ¢ = 0 denotes time immediately after load 
application 
i = value of ¢ at which lateral deflection of column be- 


comes infinite 


(422 a ) © tea ae 
¥ = e oy ) 
‘ 2h a 


YRr/V, initial-straightness parameter, e~ (4/4) l¢, FRI 
bo = built-in lateral deflection at mid-height of column 
5(1) = lateral deflection in excess of 69 at time ¢ at mid- 
height of column 
= compressive strain 
€,(t) = compressive strain at time / in left or concave flange 
at mid-height of column 
ent) = compressive strain at time ¢ in right or convex flange 


at mid-height of column 


$59 
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t 
Fic. 1. Approximation of a smooth stress history by a series of 
infinitesimal steps 
nytt) = parameter proportional to plastic strain at time ¢ in 


left flange at mid-height of column, 


EB Fi ( su) 
- €7 ws 
SK «,—8 E 


nplt = parameter proportional to plastic strain at time ¢ in 
right flange at mid-height of column, 


EB a j ‘*) 
_ ER ‘ 
2K ¢, —@& ( E 


o = compressive stress 

o1(t) = compressive stress at time ¢ in left flange at mid 
height of column 

o,(t) = compressive stress at time ¢ in right flange at mid 
height of column 

a = average applied stress on column 

o, = Euler stress of column, 7? /4(L/b)? 


CREEP EQUATION FOR MATERIAL UNDER VARYING 
STRESS 


Because the stress in any fiber of a creeping column 
varies with time, the description of the material given 
by Eq. (1) is not sufficient for analyzing a column made 
of this material. A more general description of the 
material, which includes the behavior under varying 
stress, will therefore be formulated with the aid of 
Shanley’s engineering hypotheses regarding creep. * 

Shanley’s hypotheses, for a material at constant 
temperature, may be stated as follows: 

Iypothesis (1). 
dition of stress o and strain ¢ and the stress is held con 
stant; then the strain rate (O€/O/), constant 18 a function 
only of the stress and the strain and not of the path 


Suppose the material is at some con 


whereby these conditions were attained. In conse- 
quence of this hypothesis, the strain change de occurring 
during an infinitesimal time interval d/ at constant stress 
may be written 


de = fle, €) di (2) 


ITypothests (2). 
dition of stress o and strain ¢ and the stress is given an 


Suppose the material is at some con 


instantaneous change; then the material behaves elas 
tically, with the elastic modulus / being that associated 
with the temperature of the material. In consequence 
of this hypothesis, the strain change de occurring during 
an instantaneous stress change do may be written 
de = (1/E) da 2 

A general expression for strain rate under variable 
uniaxial stress can now be formulated by considering 
the stress-time curve to be made up of a series of infin; 
tesimal steps (see Fig. 1). Each step consists of a con 
stant-stress period of time dt followed by an instan 
taneous stress increase of an amount do. The strain 
changes occurring during each of these two portions oj 
the step are given by Eqs. (2) and (3), respectively, 
and the total strain change is 


de = f(a, €) dt + (1/E) da 


The time interval consumed by the step is dt. Divi 
sion by d/ therefore gives the average strain rate, or 


de/dt = f(o, €) + (1/E) (da/dt) (4 


Eq. (4) constitutes a general creep law for a material 
under varying stress. Its application requires a knowl- 
edge of the creep function /, which may be defined as 
the strain rate under constant stress expressed in terms 
of the stress and the strain. For the material under 
consideration, this function can be obtained from Eq 
(1). Differentiation of Eq. (1) with respect to time 
gives the strain rate under constant stress as a function 
of stress and time 


1 


* ,Bao;kh t - 
(O€/Of )g constant 4 1 Ker’t v 


The time, however, is related to the stress and strain 


through Eq. (1) that is, 


Bo;1 b 
Ae”’! {} 


{= ,1€ = oe EK) 


Use of Eq. (6) to eliminate ¢ in Eq. (5) gives the strait 
rate under constant stress in terms of stress and strain 


(>) K(Ae*’)'/4 
ol [fe — (e/E)“ —' 


o constant 
Therefore, 


Kise)" 

ee on te E)|“/ 

Eq. (7) gives a formula for /(¢, €) which is particu 
larly suitable for compressive stress (@ positive). For 
tensile stress (o negative) the formula gives small, but 
still positive, values for /(o, €)-that is, the creep under 
tensile stress is still compressive in character, according 
to Eq. (7). Rigorously, therefore, the column analysis 
that follows should be stopped when the total lateral 
deflection at the mid-height is equal to half the column 
width, for at that time the stress in the convex flange 
passes from compression to tension, and a new formula 
for f(o, €) should be used for the tension flange in place 
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of Eq. (7). However, for simplicity, such a substitution 
will not be made, and Eq. (7) will be assumed to hold 
thoughout the analysis. The errors thus incurred 
are described in the ‘‘Discussion” section, where it is 
shown that they would generally be small in practical 


cases. 
ANALYSIS OF COLUMN BEHAVIOR 


Description of Problem 


The column before application of load is shown in 
Fig. 2a. It consists of two thin flanges, each of cross- 
sectional area .S/2, whose centerlines are kept a distance 
b apart by a web of negligible area but infinite shear 
stiffness. A slight initial crookedness is present in the 
shape of half a sine wave, the magnitude of the mid- 
height lateral deflection being 6). The law of material 
behavior is taken to be Eq. (4), with /(o, €) defined by 
Eq. (7). 

A load is instantaneously applied producing an aver 
age stress 6 on every cross section and an instantaneous 
static deformation as shown in Fig. 2b. By hypothesis 
the material behaves elastically during instantaneous 
loading. Therefore, the conditions existing immedi- 
ately after load application (¢ = 0) are obtainable from 
a conventional elastic analysis and, assuming small de- 
flection, may be shown to be as follows: The deflec- 
tion shape remains half a sine wave, but the mid- 
height lateral deflection is increased above 6) by an 
amount 6,(0), where 


6:(0) = d9/%/(o, — &) (S) 


and o, = m°k/4(L/b)*, the Euler buckling stress for 
an elastic column. The compressive stresses developed 
in the left (concave) and right (convex) flange at the 


mid-height of the column are, respectively, 


by 0, 
a(0) =a} 1+ 2 - (9) 


b \o, — & 
bu 0, 

or(O) = «| 1 — 2 ( (10) 
b \«a, — & 

Their difference, for future reference, is 

46, 0G, 

o1(O) — or(O) = & (11) 

b o, - @¢ 


rhe compressive strains in the left and right flange at 
the column mid-height are, respectively, 


(0) = ao (O)/E (12) 
er(O) = or(O) Ek (133) 


If the load on the column is now held constant, the 
column will continue to deform because of creep and to 
deflect laterally because of more rapid creep in the more 
highly stressed flange (see Fig. 2c). The problem is 
to determine the time history of the column deforma- 


tions. 
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(a) (b) (c) 
Fic. 2. Schematic description of column. (a) Before load 
ing. (b) Immediately after application of load (¢ = 0 c) At 


time ¢ after application of load 


Development of Basic Equations 


Requirements of equilibrium, compatibility, and ma 
terial behavior must be satisfied at all sections of the 
columns at every time ¢. For simplicity, however, the 
lateral deflection shape is assumed to be merely a magni- 
fication of the shape at ¢ = 0—that is, half a sine wave 
and the deflections are assumed small compared to the 
length. Because of this reduction to a single degree of 
freedom in the deflection shape, the requirements of 
equilibrium, etc., need be satisfied only at one section. 
Selecting the mid-height section, one may express the 
requirements in equation form as follows: For equilib- 
rium, 


ao; + or = GE (14) 
on — or = (46 6) (69 + 6) (15 


The assumptions that the web has infinite shear stiffness 
and that the deflection shape is half a sine wave of 
small amplitude imply the following condition of com 
patibility : 

(e, — €r)/b = increase in curvature at center of 


column 
6,\(3° L?) (16) 


Finally, the law of material behavior must be satisfied 
in each flange, or, from Eq. (4), 


le, — (or Bk) = /(aoL, € (17) 


(d dt) 


(d/dt) ma ie — ke) l(or, €r) (IS) 


where /(¢, €) is defined by Eq. (7 

Eqs. (14) to (18) constitute a system of five equations 
in five unknown functions of t—that is, o,, or, €1, €p, 
and 6,. These equations are to be solved subject to 
the initial conditions existing immediately after the 
load is applied (¢ = 0), given by Eqs. (8)-(10), (12), 


and (13). 


Reduction of Basic Equations 
The system of Eqs. (14) to (1S) can be reduced to 


two equations in terms of the plastic strains ¢ 
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(o.,/E) and ee — (or/F). Eq. (16) is first used to vr [{7™ are “2p vs 
eliminate 6, in Eq. (15). Eqs. (15) and (14) are then v1. J0 ii sehen J, e. =e dy (29 , 


solved simultaneously for o, and op in terms of e, — 


(o),/E) and er — (or/E) to obtain In Pearson's ‘Tables of the Incomplete l’-Function,”s 


values are tabulated of the quantity J(u, p), defined 


Kk 5 ; 
o, = o(O) + = z ~ | («: a “*) = («« - m) by 
26, —6 : ; 


| ele as 
(19) I(u, p) = | e "v’ dy 
M(p + 1) Jo 


Eq. (29) can therefore be put in the alternate form al 


(20) YR = l = | 
(m V/4K,—-— 1] = 1( aw. -,—- 1) 30 vi 
YI K K r 


Using Eqs. (19) and (20) to eliminate o, and ox in Eqs. 


or = or(O) — B- 








(17) and (18) then gives With nx known as a function of 9, through Eq. (29) or 7 
dy, /dt = Kyz,e— 8/9," —' (21) (30), Eq. (21) can be integrated to obtain ¢ as a fune 
tion of n,. The result is fi 
dynr/dt = Kypre™®—™ or si (22 7 ian 
esi 
where Ww kK | =< my 2M (31 si 
1, = oe . - (« = “) (22 where the integration can be performed numerically, al 
ee ae : Examination of Eqs. (50) and (81) provides some in ce 
KB 6 pe sight into the behavior of the quantities ny. and npr, which ti 
= — i (« ate ) (24) are proportional to the plastic strains. Because J(u, p n 
2Ko,.—¢6 E ; . eee 
tends to unity as u approaches infinity and because the st 
AEB 3 1/K aes om initial-straightness parameter Yr/y. lies between one re 
= ( a oan ) e : (-)) and zero, Eq. (30) indicates that nr approaches a con in 
stant value as n; goes to infinity. Because nz is bounded, 2 
. AEB a ) a oat a Eq. (31) then implies that y,/ is also bounded as 7, ap 
ie ( Maat . : (<0) proaches infinity — that is, in a finite time infinite strain 
is developed in the left flange. 

The problem is now reduced to that of solving the Special Forms of Solution..-For several special cases " 
differential Eqs. (21) and (22) for n, and ne subject to of interest, the solution, Eqs. (29) and (30), can be re ‘ 
the initial conditions 7,(0) = nr(O) = 0. duced further. For the extremely crooked column 

that is, for yr vy, > 0 —these equations reduce to 
Solution of Reduced Equations 
ne = O (32 

General Solution. —Eqs. (21) and (22) can be solved 
simultaneously to obtain an equation relating n, and l “my . via st 
nr and one relating ny, and ¢. First, df is eliminated WwW = K | cn mt dn 
between Eqs. (21) and (22), and the variables are sepa- F 
rated to obtain | | | la 

( a(n ee \) (59 
VR (o. \Q/K)—1 2-29, 15 ' 1/K)—1 0 —2np 17 “ 7 . 
(2n,) e ~"- d(2n,) = (2nr) e "* d(2nr) 
: For the perfectly straight column that is, for yr /7 
(27) 1 the equations become 
where the parameter yr y, may be expressed as follows os (34 
[see Eqs. (25), (26), and (11)]: : 
vit = m1" (35 
YR =e (#/KIley0) —eRO) _ 9~(8/K)(Ah/s) aee/se—3)1 (gy 
YI Ihese results correspond to direct compression, and 
Eq. (35) is equivalent to Eq. (1). 
and may be regarded as a measure of initial straight- For the special case A = 1—-that is, for a viscous to 
ness, the value unity corresponding toa perfectly straight material with a rate of flow proportional to e”’ |see Eq ar 
column and the value zero to an infinitely crooked col-  (1)} Eqs. (29) and (31) reduce to y 
umn. Integration of Eq. (27) and imposition of the Vi 
initial conditions yield the following relationship be- —— | he [ 7; | (36 as 
tween n,(¢) and nr(t): ™ » YI pt 
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Graphical Representation of Solution. Plots of 1; 
and nr against the time parameter y,/ are given in Fig. 3 
for four values of the material constant A and several 
values of the initial-straightness parameter yr, Y. 
The plots confirm the qualitative deduction that, for 
a given value of yr/yv less than unity, nr approaches a 
constant as 7, approaches infinity. Except for the 
perfectly straight column (yr/y, = |), n, becomes in- 


finite in finite time. 


Stresses, Strains, and Lateral Deflections 

The quantities n, and nx, defined by Eq. (29) or (30) 
and Eq. (31) and plotted in Fig. 4, are the key to the 
computation of the stresses, strains, and lateral deflec- 
tions at the mid-height of the column. ‘The quantites 
n, and nr are themselves proportional to the plastic 
strains in the two flanges, as Eqs. (23) and (24) indi- 
cate. In terms of n, and ne the other quantities of 
interest are expressible as follows: From Eqs. (19) and 


20), the flange stresses are 
ao, = 01(0) + (AB) (ni — npr) (38) 
or(O) (KB) (n nr) (39) 


Or 


From Eqs. (23) and (38) and (24) and (39), the strains 


a,(Q) AK (=" ’) 
, (40) 
ee sc Se OF Aoi 
or(0) ' K (™*: ’) ] ‘i 
= 2" oe "OE ies Mia 


Finally, from Eqs. (16), (40), and (41), the ratio of the 
lateral deflection 6, to the column width d is given by 


6 1 | o,(0) — opx(0) ' K 
= (n, - ) 
b 4 ‘ —- (= 


=| 5 ) _& 
— fon <= . 
i\e <2 —_— 7 


t 


are 


6,(0) K 

b * 2B ™ _ 

As Eq. (42) indicates, 9), — ne is directly proportional 
to the lateral deflection due to creep. In Fig. 4 plots 
are therefore given to show the variation of nL — nr with 
yit for four values of K and several values of Yr/YL- 
Vertical dashed lines are placed in Fig. 4 to indicate 
asymptotic values of y,t as the lateral deflection ap- 


proaches infinity. The value of ¢ corresponding to an 
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Fic. 3. Variation of plastic-strain paramaters n, and n, with 


time parameter y,/ for various values of the initial-straightness 


parameter y,/¥, 


asymptote represents an upper limit to the lifetime of 
the column and will be denoted by ¢,,. 

In most cases the column lifetime will be of primary 
interest. In Fig. 5, therefore, the asymptotic values of 
yit, denoted by yz/,,, are plotted against the initial 
straightness parameter yz, y, for the four values of K 
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Fic. 5. Variation of lifetime parameter y,¢.. with initial 
straightness parameter y,/y, for four values of the material 
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stant A for very crooked columns (Y¥p/¥, = Y) is 
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considered in the preceding graphs. The limiting vari- 


ition of y./-- with A for an extremely crooked column 
(yey. > 0) is shown in Fig. 6. 


AppLICATION OF RESULTS TO SPECIFIC MATERIALS 


The foregoing general results may be applied to a col- 
umn made of a specific material if the constant-stress 
creep curves of the material have the form defined by 
Eq. (1) and if the constants in this equation are known. 
The column behavior thus deduced will be accurate to 
the extent that the other assumptions are correct. 


Data on Material Constants 


Some information has been published which indicates 
that Eq. (1) may be representative of a number of ma- 
terials. Reference | suggests the use of Eq. (1) for 
75S-T6 aluminum alloy at 600°F., stabilized at this 
temperature for 100 hours prior to use, and gives the 


following values for the constants: 


E = 5.2 X 10° Ibs. per sq.in. 

K = 0.66 

A = 2.64 X 107 (¢ measured in hours) 

B = 1.92 X 10~* (o measured in Ibs. per sq.in.) 


Fig. | of reference 2 indicates that Eq. (1) also applies 
to a low-alloy steel (0.40 C, 1 Cr, 0.5 Mo, and 0.25 V) 
at 800°F., which, prior to test, is heated as 3-in. round 
for 5 hours at 1,650°F., oil-quenched, and tempered 
for 12 hours at 1,240°F. and then for 5 hours at 1,100°F. 
The constants for this material are as follows: 


E = 26.5 X 10° Ibs. per sq.in. 

K = 0.25 

A = 6.24 X 10~* (¢ measured in hours) 

B = 7.27 X 10~° (o measured in Ibs. per sq.in.) 


A limited amount of data in reference 2 suggests that 
Eq. (1), with the following values of the constants, may 
also apply to the same steel at 1,100°F.: 


k= 22.2 X 10° Ibs. per sq.in. 

K = 0.29 

A = 31.3 X 10~¢ (¢ measured in hours) 

B = 30.7 X 10~° (o measured in Ibs. per sq.in.) 


Column Curves 


The information on column lifetime is given rather 
compactly in Fig. 5. When dealing with a specific 
material, however, it may be desirable to put the life- 
time data in the more familiar form of column curves. 
Typical curves of this type are presented in Fig. 7 for 
the low-alloy steel at SOO°F. Each column curve 
applies to a particular crookedress ratio 6)/L and to a 


specific lifetime ¢, 
DISCUSSION 
Significant Aspects of Results 
Existence of a Critical Time.—The present study re- 


veals the existence of a ‘‘critical time’’—that is, a finite 
time at which the deflections approach infinity asymp- 
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Column curves for H-section columns of a low-alloy 
steel at 800° F 


totically (see Fig. 4). This behavior is a consequence 
of the exponential nature of the creep function f [see 
eq. (7)] and is not generally to be expected for mate 
rials whose creep rates are less severely affected by 
In particular, a careful study of the linearly 
that is, a column for which /(¢, €) 
by Prof. Joseph 


stress. 
visco-elastic column 
takes the form o/\ with \ constant 
Kempner and his associates at the Polytechnic Institute 
of Brooklyn,° has shown that no critical time other than 
{ = o exists for such columns. 

The critical times defined by the asymptotes in Fig. 
t and summarized in Fig. 5 have practical design signifi- 
cance provided the lateral deflections do not become 
excessively long before the critical times are reached. 
A numerical check of the lateral deflections correspond- 
ing to the steeply rising portions of the n, — nr curves 
of Fig. 4 (say, at n, — nr = 6) and to times close to the 
asymptotic value /,, indicates that, in this region and 
for columns of reasonable proportions, the lateral de- 
flections are generally still less than the column thick- 
ness. Hence, ¢,, could be expected to approximate 
closely the actual column lifetime were all the other 
assumptions of the analysis correct. 

Effect of Initial Crookedness.—That initial crooked- 
ness is theoretically essential to creep buckling is indi- 
cated in Fig. 4, which shows that n, — nr remains zero 
at all times when the column is perfectly straight 


(yr/y_ = 1). No lateral deflections develop, and the 
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column merely behaves as a compression creep speci- 
men in accordance with Eq. (1). 

The sensitivity of yzf,, to the initial-straightness 
parameter yr/yz. is shown in Fig. 5. The dependence 
of ¢,, itself on 69 is seen to be even more acute than that 
of yzter ON Yr/¥_ When one realizes that y, varies 
directly, and yr/y. inversely, to e raised to a power 
directly proportional to 4p. 

A comparison of parts (a) and (b) of Fig. 7 shows 
graphically the large effect on ¢,, of doubling 40, es- 
pecially at the lower values of L/b. 

Effect of Applied Load.—The extreme sensitivity of 
column lifetime /,, to the applied load can be revealed 
by a study of the parameters y, and yr/y. and Fig. 5. 
However, Fig. 7 is much more effective in showing this 
sensitivity. Either part of Fig. 7 shows that only a 
slight change in @ is needed to change /,, by a factor 10. 


Extension of Results 


Design situations may involve eccentricity of loading 
é) combined with an initial crookedness 69. If this case 
is analyzed in a manner similar to that described for the 
case of initial crookedness alone and the assumption is 
made that successive deflection shapes are all affinely 
related to the deflection shape existing immediately 
after application of load, then the results are identical 
with those already obtained, provided that in the latter 
dy is replaced by 6) + e and go, is interpreted as the 
Euler load corresponding to the deflection shape exist- 
ing after load application. It will always be conserv- 
ative, however, to use for o, the value corresponding to a 


sinusoidal shape. 


Validity of Assumptions 


Material Behavior.—Few if any general truths are 
available regarding the creep behavior of materials 
under varying stress. Shanley’s theory* was selected 
as the basis of the present analysis because of its sim- 
plicity and in the hope that, despite its admitted short- 
comings, it could yield results sufficiently accurate for 
some engineering purposes. 

The most questionable aspect of the theory is prob- 
ably the assumption that the creep rate under constant 
stress is a function only of the instantaneous values of 
the stress and strain and not of their histories. The 
assumption is questionable in view of some scattered 
experimental evidence that indicates that history does 
affect the constant-stress creep rate of alloys. The 
magnitude of this effect has not been systematically in- 
vestigated. Where the effect is large, the present anal- 
ysis would not apply. 

The assumption of elastic behavior under instantane- 
ous application of load is better founded than the pre- 
viously discussed assumption, particularly at higher 
temperatures. At sufficiently high temperatures, ex- 
periments indicate the nonexistence of a unique stress- 
strain curve and no clear demarcation between plastic 


flow and creep. The stress-strain curves seem to ap. 
proach the elastic line as the rate of load application 
goes up. Further evidence to support the assumption 
of elastic behavior under rapid loading has been ob. 
tained from tests in which creep at constant stress has 
been followed by a sudden small change in stress 
In such tests, the ratio of the sudden stress change to 
the resulting strain change is observed to be approxi- 
mately equal to the initial elastic modulus. The 
phenomenon of elastic behavior under rapid loading 
could explain the results of experiments intended to 
confirm the theory of propagation of plastic waves jn 
bars.6 These tests involved the sudden addition oj 
an increment of tensile stress in a bar already stressed 
in the plastic range. The wave-front velocity corre- 
sponded to the initial elastic modulus rather than to 
the tangent modulus. 

Apart from questions of the basic creep theory em 
ployed, some attention should be given to the conse- 
quences of the particular formula { Eq. (7)] taken to de 
fine the constant-stress creep rate function f(¢, €). This 
formula implies that the creep under tensile stress 
(negative values of o) remains compressive (positive 
in nature. Because of the nature of e raised to a neg- 
ative power, however, the amounts of compressive 
creep given by the formula will be small enough to be 
considered zero, and it may be stated, therefore, that 
the present analysis is effectively based on the assump 
tion that no creep occurs in tension. The analysis is 
consequently rigorously valid only up to the point where 
the tensile stress on the convex side of the column be 
comes large enough to cause appreciable tensile creep, 
and continuance of the analysis beyond this point must 
result in an overestimate of the column lifetime. The 
correct column lifetime is therefore bracketed between 
a lower limit (the time at which the analysis becomes 
invalid—-say, conservatively, when org = OQ) and an 
upper limit (the computed column lifetime). 

Despite the neglect of tensile creep, physical con 
siderations indicate that in most practical cases the two 
limits are close together—that is, the analysis does not 
become invalid until almost all of the computed critical 
time has passed. This is a consequence of the fact that, 
when the analysis becomes invalid (that is, when or = 
0), the stress o, in the concave flaige is twice as high as 
the mean column stress ¢. Therefore, if & is sufficiently 
high to imply efficient use of the material [but not so 
close to the Euler stress that, in conjunction with the 
initial crookedness, it causes the initial stress ap(0) to 
be far below @], then by the time the analysis becomes 
invalid o,; will be in the region of extremely rapid creep 
and collapse of the column will be close. Consequently, 
the computed critical time may be taken as the correct 
critical time with little error. 

Numerical checks confirm the contention that for 
practical cases the time at which or = 0 and the com- 
puted time of collapse are close together. In cases of 
doubt, both of these times should be computed, and, 
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where they differ greatly, judgment should be used in 
lecting some intermediate value as the correct column 
In many cases a lateral deflection equal to 


se 
lifetime. 
one-half the column width would be considered exces- 
sive and taken as failure; in such cases, of course, the 
present analysis is valid over the entire time range of 
interest. 

Column Shape.—The assumption that the deflection 
shape remains at all times affinely related to the sinu- 
soidal shape existing immediately after application of 
load is known to be incorrect. The deflection shape 
may be sinusoidal immediately after load application, 
but, as tests of rectangular-section columns indicate, ! 
near failure the column has a high curvature (nearly a 
kink) in a short central portion and is relatively straight 
elsewhere. However, the assumption of an unchang- 
ing column shape is felt to be justified in view of the 
uncertainty regarding material behavior and the great 
simplification it results in. 


CONCLUDING REMARKS 


On the basis of certain assumptions regarding the 
column shape and the creep behavior of the material 
under varying stress, an analysis has been made of a 
slightly crooked pin-ended H-section column carrying 
aconstant load. The applicability of the analysis to at 
least two materials (75S-T6 aluminum alloy at 600°F. 
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and a low-alloy steel at 800° and 1,100°F.) has been 
indicated. 

The.analysis results in information regarding the 
column lifetime and shows the sensitivity of column 
lifetime to initial crookedness and load. The results 
can be put in the form of column curves giving, essen- 
tially, strength (for a given initial crookedness and de- 
sired lifetime) vs. slenderness ratio. 
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Buckling of Sandwich Cylinders Under 
Bending and Combined Bending and Axial 
Compression’ 


CHI-TEH WANG? anv D. P. SULLIVAN? 
New York Unwersity 


SUMMARY 


A theoretical investigation of the buckling of sandwich cylin 


ders under bending and combined benditig and compression is 


carried out. 


the effects of transverse shear using Galerkin’s method, the buck- 
For a sandwich cylinder with a weak 


ling load is determined. 
core, it is found that buckling occurs when the maximum com 


pressive stress in the cylinder reaches the buckling stress of the 


cylinder under axial compression alone 


SYMBOLS AND UNITrs 


a = radius of cylinder to mid-plane, in. 
c = transverse shear stiffness of isotropic sandwich 
cylinder, C = hG, lbs. per in. 
D = flexural stiffness of isotropic sandwich cylinder, 
in.lbs. D = Eth?/2(1 — v?) 
E = Young’s modulus for the faces, lbs. per sq. in. 
G = shear modulus, G = £/2(1 + v) 
h = depth of isotropic sandwich plate measured between 
middle planes of faces, in. 
I = moment of inertia of the sandwich cylinder about 
its diameter 
l = length of cylinder, in. 
M = bending moment, in.Ibs. 
m = number of half waves longitudinally 
N,, N, = resultant normal forces in x- and y-directions, in. 
Ibs. per in. 
1, ae = resultant shearing force in xy-plane, lbs, per in. 
N, = applied normal force in axial direction 
n = number of full waves circumferentially 
t = thickness of the faces, in. 
u,v, U displacements in x-, y-, 2-directions, respectively, 
of a point in middle surface of cylinder, in 
x,¥,2 = rectangular coordinates (Fig. 1) 
v = Poisson’s ratio for the face material 
7] = circumferential coordinate 
( o? 2 ) 
did = operator, i ‘ - 
Ox? Oy? 
(= oe mae “)-(2 : ) 
ox? ~ Ox? oy? oy? Ox? a?*0oe? 
( °? o? ) (= o? ) 
a = operator, a ,| = isp" meer 
Ox* oy Ox* a*og- 
A = mra/l 
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By solving Donnell’s equation modified to include 


INTRODUCTION 


N PREVIOUS PAPERS,' + the buckling of sandwich 
FF ciate under axial compression and_ torsional 
loads have been investigated. In the present paper, 
the investigation is extended to the buckling of sand 
wich cylinders under bending and combined bending 
and axial compression. By solving Donnell’s equation 
modified to include the effects of transverse shear for 
sandwich curved plates and shells following the Galer 
In the 
case of sandwich cylinders with weak core, the buckling 


kin’s method, the buckling load is calculated. 


load is found to be only a function of the core rigidity, 
which agrees with the test results in the limiting case 
of sandwich cylinders under axial compression. ! 


FORMULATION OF THE PROBLEM 

By assuming isotropic core material and neglecting 
the bending rigidity of the faces about their own middle 
surface, the equilibrium equations for an element of a 
sandwich cylinder have been derived in reference 2. 
By neglecting terms that were regarded as small by 
Donnell,® these equations can be reduced to a single 
equation in terms of the lateral deflection w only 
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This is the so-called Donnell’s equation modified to in- 
clude the effects of transverse shear for sandwich curved 
plates and shells first obtained by Stein and Mayers.°® 

In the case of sandwich cylinders under combined 


bending and axial compression, we have 


N, = Ny = 0 
N, = —N, + (2Mat/I) cos 6 (2) 


where .V; is the force per unit length due to axial com- 
pression, M is the bending moment, and / is the mo- 
ment of inertia of the cylinder. With such a loading, 


Eq. (1) thus becomes 


D 
DViw + (: — 6 v°) 


j2tk O*w , : 2Mat O-w || 
. _ y* —N, + cos 6 : ( = () 
(a? Ox! I Ox? 


(:3) 


Eq. (3) may be solved by Galerkin’s method as 
follows: First, the deflection w is assumed in the form 
of a series that satisfies the boundary conditions but with 


undetermined parameters. Assume that the cylinder 
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is simply supported at the ends. The boundary con 


ditions are then w = Oand \/, = Oatx = Oandx = /. 
With terms neglected by Donnell, these conditions be- 
come w = Oand 0°w Ox? = Oatx = Oandx = It is 


easy to see that the series 


MTX C 
, me 


, cos né (4 


w= sin 


satisfies these conditions where C, are undetermined 
parameters. 

Next, substitute the assumed series for w, expression 
(4), into Eq. (3). If expression (4) happens to be the 
exact solution of Eq. (3), after substitution Eq. (3) will 
be identically equal to zero. In general, this is not so, 
and the resulting expression will be a function of x and 
6 which we shall denote by Q. Galerkin’s equations for 
the determination of the coefficients C,, are 


_ mmx ‘ 
Q sin cos n6 ad@ dx = 0 >) 


,1,2,...). Substituting, Eq. (5) becomes 


O-w D ~) 
‘ x 


> O'w DF_.dOtw\ . mrx ’ Fe 
DV' w+ F =~ —— Y= ) sin cos n6 ad@ dx + Ny Vv —-— Ye 
0 Jo ( Oox4 ] J 0 Ox? Ox? 
_ mrx 2Mat [{* [' 
sin i cos n6 ad6 dx — F 


O-w D , O-w Manx ; 
— Vilcos 0 ox? + C V®l cos 6 > 4 x sin —— cos 6cos nO ad6 dx = 0 (6) 
: od 4 l 


Ox 4 


x*F J 


where F = 2Et/a®. Carrying out the integration and introducing the notation \ = mma /, Eq. (6) becomes, after 


dropping a common factor, 


_J 2 2)4 i 1,4 FD 2,4()2 4 2 Ni) aye 2\21 [a2 4 D v 2 
Cn (A vey + e8 +3,7°"¢ Te) - > (A* “> 98°) Ft 7 ¥ nN { 


Mat Jf .. 7 a D F eee. ; 
Mi (C, 1 (a: _ C 6n } (A* + n*)* + C [(A? + n?)? 


ID 


2n(6A\2 + 10n? + 3) + (3A? + 15n? 4+ 1)) 


Po : i eae i 
C41) ato 62) (A? + nw) + GIO + m*)* +- SA 


(3A2 + 15m? + 1)] 


form ~ 1 where C_,; = 0. 


+ qg27{/(2n? + 6n? + 1) — 4n(A* + n° 4+ 21) ea 


+ q@?/(2\° + Gn? + 1) + 


{ 


+ 


+ 3(A2 + n*)? + 127(\? + n*) — 


24 2)? + 12m2%(A2 + un?) + 2n(6dA2 + 10 mn? + 3) + 


tn(r\° + n- + | 


For n = 1, 2C) should be substituted for C. 


The substitution into Eq. (6) of n from 0 to © produces an infinite number of algebraic simultaneous equations 


Since there is no constant term in these equations, one solution of these equations is obviously the trivial one 


namely, 


To find the nontrivial solution, we must set the determinant of the coefficients of the parameters C,, in these equa 


tions to zero, which gives us the buckling load. 
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BUCKLING UNDER PURE BENDING 


Let us first study the case of pure bending. In this 
case, we shall let V; = Oin Eq. (7). It is noted that, in 
Eq. (7) for a given n, C,4; always appears as an un- 
known. As a result, for any 1, we have m equations 
with (7 + 1) unknowns, and consequently these equa- 
tions cannot be solved. To avoid this difficulty Fliigge’ 
has proposed to drop this extra unknown in the last 
equation and then to calculate the buckling load by 
setting the determinant of the coefficients of C,, to zero. 
The error introduced by this omitted term becomes 
smaller and smaller as more and more parameters C,, 
are taken. The justification of the method can also 
be seen by observing that, as m increases, the buckling 
load converges rapidly. 

To study the effect of \ on the buckling load, we shall 
follow Fliigge’s method and calculate the buckling 
load of a given sandwich cylinder for various values of 
n. Since we are interested in sandwich cylinders with 
weak cores, this is calculated for such a sandwich cylin- 
der and plotted in Fig. 2 with m = 1, 2, 3,4, 5. It was 
found that, for all these values of m, minimum J/ occurs 
at\ = «. Thisis an agreement with the experimental 
and theoretical results of reference | for the limiting case 
of axial compression. The convergence of /D/2Mat 
as m increases can easily be seen from the figure. 

Having established the fact that, at buckling, \ must 
be extremely large for sandwich cylinders with weak 
cores, a simpler procedure can be followed to find the 
buckling load. 

By assuming that A is a large number, many terms in 


Eq. (7) can be neglected, and Eq. (7) becomes 

C, + (Mai/CI) (C,~1 + Cysi) = 0 (S) 
Divide Eq. (8) through by the factor (—J\Jat/2C/), 
and letting 


L = CI/Mat 


Eq. (8) may be written in the form 


Cn ss £6. + Os = © (9) 
forn~1. Forn 0, we have 
LCo T Ci = (0 (10) 
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and, forn = 1, we have 
2C5 + LC, + Co = QJ (11 


To calculate the buckling load, instead of solving the 
determinantal equation where we must drop the term 
C41 in the last equation, a more rigorous solution can 
be obtained as follows: We observe that Eq. (9) is a 
finite difference equation in terms of C,,. 
solution of Eqs. (9), (10), and (11) is 


The trivia] 


Cc =0 iw=e3681.23... 
The nontrivial solution of Eq. (9) is 
C, = Aa" + Ba 


se = 2, 2, By. (12 


where A and B are arbitrary constants and a is related 
to L by the formula 

a+aL+1=0 (13 
or 

L= —(a+ a) (14 


Substituting Eq. (12) into Eq. (11) and using relation 
(13), we find 


2C) = —LC “a C2 = —L (Aa + Ba 1) — 
B(Aa®? + Ba?) = —A (a? + al) —- 
B(a-*+a"'L) =A+B (15 


Now, Eq. (10) is the only one yet to be satisfied. Sub- 
stituting Cy given by Eq. (15) and C; given by Eq. (12 
into Eq. (10) and using relation (14), we obtain 


(A — B) (a — a") =0 
from which we find a? = 1 and a = +1. Hence, 
I, = £2 or 
2Mat/I = +C = +hG (16 
We note that, when a = +1, the two roots of Eq. 


(13) are equal and Eq. (12) must assume the form 
C, = (A + nB)a’ 
But C,, is finite when m — o: we have, therefore, B 
= (. In this case, we find that a = +1 is again 
the correct solution. The plus or minus signs in Eq. 
(16) indicate that a positive or negative bending mo- 
ment on the cylinder will make no difference so far 
as buckling is concerned. This is equal to the buck- 
ling load when the cylinder is under axial compression 
only. Eq. (16) indicates that the sandwich cylinder 
buckles when the maximum bending stress reaches the 
buckling stress when the cylinder is under axial com- 


pression only. 


BUCKLING UNDER COMBINED BENDING AND 
COMPRESSION 


The determination of the buckling load for a sand- 
wich cylinder under this combined loading can be car- 
ried out ina similar manner. Again, by assuming large 
\, Eq. (7) becomes 

(Concluded on page 485) 
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The Effect of Nonuniform Surface 
Temperature on the Transient Aerodynamic 
Heating of Thin-Skinned Bodies 


A. E. BRYSON* anp R. H. EDWARDS* 
Hughes Aircraft Company 


ABSTRACT 


A quasi-steady analysis is given of the transient forced-con- 
vection heating of bodies with thin metallic skins, where the 
boundary layer is laminar, taking into account the effect of the 
nonuniform surface-temperature distribution on the rate of heat 
transfer. Lighthill’s approximate theory of heat transfer to a 
nonisothermal surface is used and leads to a formulation of the 
problem in terms of a partial integrodifferential equation that is 
solved in closed form. The results are compared with previous 
analyses in which the theoretical heat-transfer coefficient for an 
isothermal surface is used, and it is shown that this latter approxi- 
mation leads to results in fairly good agreement with the more 
exact solution. However, if heat-transfer coefficients are calcu- 
lated from the more exact solution, they may differ greatly from 
the ones predicted by the isothermal surface theory. This indi- 
cates that caution should be exercised in interpreting experi- 
mental data obtained by the transient method from a noniso- 


thermal surface where the boundary layer is laminar. 


INTRODUCTION 


— LOCAL RATE OF HEAT TRANSFER through a 
laminar boundary layer to a semi-infinite flat 
plate with a uniform surface temperature is well known 
from boundary-layer theory. For the case of u and k 
proportional to the temperature and constant C,, it 
isindependent of Mach Number and is given by 


q = 0.332(U/vx)'*f(o)k(T. — T,) (1) 
= h(T, — T,) 
where 

gq = heat flowing through a unit surface area in a 
unit time 

l’ = fluid velocity just outside the boundary layer 

v = kinematic viscosity of the fluid 

h = heat-transfer coefficient 

x = distance from leading edge of the plate 

k = thermal conductivity of the fluid 

Tl. = T. + r(a) (U?/2C,) 

Tl, = fluid temperature just outside the boundary 
layer 

C, = specific heat per unit mass at constant pressure 
of the fluid 

I’, = surface temperature of the plate 

¢ = C,u/K = Prandtl Number of the fluid 

# = viscosity coefficient of the fluid 


Presented at the Aerodynamics Session, Twentieth Annual 
Meeting, I.A.S., New York, January 28-31, February 1, 1952. 
* Research Physicist, Missile Aerodynamics Department. 
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and f(a) and r(o) are well approximated by o ° and 
a’, respectively. 

Chapman and Rubesin' have given a solution to the 
boundary-layer equations for heat transfer to a flat 
plate with a surface temperature that is expressible as 
a power series in the distance from the leading edge 
Lighthill’? has subsequently presented an approximate 
solution for heat transfer to a body in terms of the 
skin friction and temperature distributions along its 
surface; in Lighthill’s solution the temperature dis- 
tribution need not be analytic as in the solution of 
Chapman and Rubesin. For the flat plate, Lighthill's 
solution (with the constant multiplier reduced 2 per 


cent) is 


: “* d[T,(&) — T,] 
gq = —0.332(U/vx) “* ak : 2) 
Fs oll — (&/x) “* 


$ 


It is easily seen that Eq. (2) reduces to Eq. (1) for the 


case 7, = constant, since the integral appearing in 
Eq. (2) is a Stieltjes integral, and 7, — 7, = Oforé = 
0. If 7, is not constant but the denominator of the 


integral of Eq. (2) is approximated by unity, then 
Eq. (2) also reduces to Eq. (1); this is essentially the 
approximation that is made when one uses the isother- 
for surfaces that are 


mal ‘‘heat-transfer coefficient”’ 


not isothermal. 


Although both Eqs. (1) and (2) are strictly appli- 
cable only to steady-flow problems, it is possible to 
utilize them as in unsteady-flow 
analyses with certain 
indicated that, if the quantities (x"/U"*') (d"*l/dt' 


, are small compared to unity, the 


approximations 


restrictions. Moore has 


n = il, 2, 
boundary layer on a body at any instant is well 
described by the steady-flow results that correspond to 
the instantaneous velocity and fluid properties—i.e., 
the boundary layer may be considered to be quasi 
steady. Evaluation of these quantities for most pres- 
ent-day aircraft and missiles indicates that they are 
indeed small for most technically significant condi- 
tions. 

Eqs. (1) and (2) may also be applied to a wedge if 
the bow shock wave is attached to the leading edge 
The fluid properties that are used in this case are those 
just outside the boundary layer. 
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(GENERAL PROBLEM 


The general problem with which we are concerned 
is that of predicting the temperatures of a body that 
moves with variable velocity through an atmosphere 
whose properties vary with altitude. The temperature 
of the surface will then be a function of time. Con- 
tinuity conditions require the local heat flow rate from 
the boundary layer to the surface to be equal to the 
local heat flow rate through the surface of the body 
plus the net radiation from the surface. This latter 


heat flow rate is given by 


—k,(O7/On) — (radiation to surface) + 


q = 
(radiation from surface) (3) 

where k,, = thermal conductivity of the body material 

and » = inwardly drawn normal distance from the 


surface of the body. 

Before Lighthill’s analysis was available, it was com- 
mon practice to approximate the local heat flow, q, by 
Eq. (1), despite the fact that the surface temperature 
was nonuniform. The problem was then reduced to a 
nonsteady heat-conduction problem in the body, with 
the boundary condition obtained by equating the right- 
hand sides of Eqs. (1) and (3). A more exact solution 
would be obtained by using Eq. (2) rather than Eq. 
(1) in defining the convective boundary condition. 
Such a formulation would take into account the effect 
of the variation of surface temperature on the heat 
transfer, but the problem would usually be so compli- 
cated that numerical methods of solution would be 


necessary. 


APPROXIMATE FORMULATION FOR A BODY WITH THIN 
SKIN 

For the special case of a body with a thin metallic 
skin, the heat flow normal to the surface is sufficiently 
large that the temperature is effectively constant 
through the thickness. Thus, the heat flow into a unit 
surface area of the skin neglecting radiation, can be 
approximated by 


q = G(OT, Ot) (4) 


where G = specific heat per unit surface area of the 
skin (G = Cp,,6, where C = specific heat per unit mass 
of the skin material, p,, = mass density of the skin ma- 
terial, and 6 = skin thickness).* If the right-hand 
sides of Eqs. (1) and (4) are equated, there results a 
simple ordinary differential equation for 7’, as a func- 
The distance from the leading edge, x, 
Lo? has previously formu- 


tion of time. 
enters only as a parameter. 
lated and solved this problem for arbitrary time vari- 
ation of h. A more exact solution would be obtained 
by equating the two values of qg given by Eqs. (4) and 
(2). This operation yields the following integropartial 

* This neglects the effect of heat conduction along the skin 
which 1s not completely justifiable in some cases; see Appendix 


(1) 


THE AERONAUTICAL 


SCIENCES UU & %, 


differential equation: 


G Ob 9 _ Alt) £* d[T.(é t) — T(t) 
Or ae | o [1 — (&/x)"] . 


where A(t) = 0.332(U/v)'?0'k. 


SOLUTION FOR THE STEP FUNCTION IN VELocITy 


Although this quasi-steady theory is not strictly ap- 
plicable to the case of a step function in velocity, it is 
instructive to solve this case, since it serves as a tool 
in solving the problem for arbitrary A(‘) and 7.(1), 
Since the heat-transfer rate tends to infinity as x ap- 
proaches zero, it is evident that 7, = 7, at x = 0 for 
allt > 0. We therefore consider the problem in the 


following form: 


OT, (é, 1) 


au lé 
; G OF 7 —A I. OE 
"ora? Pao ba | 
X g fi = ‘) | (6 


X 
T.(0, t) = T, 


T(x, 0) = T; 


==, 


where A, 7, and 7; are constants. Now, Eys. (6) are 
invariant under the transformation x’ 
Hence, a solution to these equations of the form 7, = 


This solution was found to be 


, 


= ox, tf = al. 


T,(t/x'*) was sought. 
[see Appendix (2) for method of solution]: 


me i l 7195/3) (r/xt/2 — - 
7 ?_— oo rig exp (—u’*)du_ (4a) 
i is = i ie rd 3) oO ' 


T,-—fT, r{ (05/3) (r/x 


= — (7b 


where 


m 
‘{m; n} = | exp (—/)"—' dt 
0 


and + = At/G. 

It is interesting to compare this solution with that 
obtained by neglecting the effect of temperature gra- 
heat transfer 

The problem 


dients—i.e., by approximating the 


through the boundary layer by Eq. (1) 
then becomes 


G(dT,/dt) = (A/x 
T,(0) = T; f 


The solution to Eq. (8) is 
T, —T; T 
= —=]-— exp (— (9a) 
T, — 7; x 


or 


= = (9b) 


The form (9b) compares with the form (7b). 
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HEATING OF THIN 

Eqs. (7) and (9) and the rates of heat transfer cor 
responding to these two equations are plotted in Fig. 1. 
Values of the incomplete gamma function were ob- 
tained from reference 5. In Fig. 2, the heat-transfer 
coefticient defined by 


h = qx/k(T, — T, 


is plotted as a ratio to the isothermal surface heat 
transfer coefficient. It will be noticed that the dis 
crepancy is fairly large here, whereas the discrepancy 
between the two theories in Fig. 1 is not so great—.e., 
the heat-transfer coefficients are quite different, but 
the more interesting quantities, the surface temper- 
atures, are not so different. 

SOLUTION FOR ARBITRARY VELOCITY VARIATION AND 
INITIAL TEMPERATURE 


For this problem, the integral equation and boundary 
conditions are identical to those of Eqs. (6), except that 
Aand 7, are functions of ¢, while 7’; is a function of x. 

Because of the linearity of the problem, the solution 
may be expressed as the sum of the two solutions that 


satisty the following boundary conditions: 


T,(0, t) = Tb) 10 
(10) 
T(x, 0) = 0 { 
T,(0,t) =O 
(11) 


It is convenient to introduce the variables 7 and 3z, 
which are defined through the equations of transforma- 


tion, 


= x (12) 
i Coe T,’(, 7); Tix) = T,') 
LA) = Toes 
Eq. (5) then takes the form 
OT ."(é, 7) 
< : oe ag , 
O7.’(z, 7) | 2 of (13) 
Or ; - L — (/s)"7"" 


This equation is essentially of the same form as that 
solved in the preceding section. It is invariant under 
a translation of tr. A particular solution is 


| exp {| —(T'(5/3)a] 
0 


By means of a 


h(s, r, a) = ‘du (l4a 
where @ is an arbitrary parameter. 
Superposition integral involving the parameter a, a 
further solution is found to be 


i se » 3 
NZ, T) = / u(r — us)exp {—[T(5/3)ul 7} du (14b) 
J0 


Furthermore, if 7 is a piecewise continuously differ- 
entiable function, it may be demonstrated that 
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lim A(z, 7) = v(7) 


z—0 


It follows that if v(7) is defined as 7,,’(7), the function 
defined in Eq. (14b) will satisfy Eq. (5) and will fulfill 
the boundary conditions of Eq. (10). Hence, the solu- 
tion corresponding to these boundary conditions is 
given by 


* 


rT/Z 
y T,/(r — uz) exp {— [P'(5/3)u] 7} du 
0 


(15) 


T,'(3, +) = 


In a similar fashion, one may show that a solution 


of Eq. (5), which has the property that 7(0, f) = 0, 


is given by 


T.'(s, 7) = | v(r — us) exp | — [I'(5/3)u] 


To satisfy the condition that 7(x, 0) = 


| du 
(16) 
T(x) for 


equivalently 7’(z,0) = 7;’(z)], we must have 


7 


T,'(z) = / v(—uz) exp {— [T(5/3)u] 7} du (17) 
J/0 


If 7',’(s) is representable in the form 


T,'(2) = > dy’ 
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where the n’s are not necessarily integers, we conside 


v(a) defined in the form 


via) = 


> Cpa” 
n 


Since 
*o j 5 ie 2T[(2n + 2)/3 
u" exp) —|T u > du = z 
4 \ 3 f 3 (15/3) 


We find that Eq. (17) is satisfied, provided 


rts. tT) = 


. 


ing J 5) 
/ T(r — uz) exp . r (°) u 
0 0 


The corresponding solution, which is obtained from 


layer heat-transfer rate, may be easily obtained. It is 


*r/z 
f T,/(r — uz) exp (—1) 
0 


The integrals appearing in the series expressions of Ex 


T,'(2, r) = 
incomplete gamma functions. 


EXAMPLE PROBLEM FOR VARIABLE VELOCITY 


Consider the case of an exponential decrease of ve- 
locity with time at constant altitude with the initial 
plate temperature constant and less than the initial 
equilibrium boundary-layer temperature. The skin 
temperature will then first increase with time and then 


decrease. Let 


U = U,;, exp (=? ty) 
T(z, 0) = 0 
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r 3 [P(5/3)}"+! 
Cn = (-—1)’ = 
2 P[(2" + 2)/3] 
Finally, then, the solution to the problem with the 
boundary conditions 


j = Bs (°) i'(5/3) |" 6,3" 
ioe =? r 
i 2 i r(2n + 2)/3) * 


/ u” exp ‘ — E (°) “| du (19 


the formulation that uses Eq. (1) to define the boundary. 
by,2' 


du + ; uw” exp (—u) du (20 
Pie = 0) Fy. 


js. (19) and (20) can be expressed in terms of complete and 


. 


where 


and 


Finally, then, 
/ ) / a or x 
J/0 2C, 
G(r — us) ]*\ . °/ ts 
24 ty ( exp ) "= (° u { du 


This can be integrated in terms of incomplete gamma 
functions. The result is plotted in Fig. 3 for the cas¢ 
Tr, = 7 


proximate Eq. 


>, along with the result obtained using the ap 
(20). 


, Gz _,. (ox an 'G 
z= = 1.905 ( o 
2A ito l hty 


Fig. + shows the heat-transfer coefficients calculated as 


The variable z’ used is defined 


by 


in the case of the step function in velocity. It is seen 


that they differ greatly from the ‘isothermal heat- 
transfer coefficient,’ despite the fact that the tem- 
peratures calculated by the two methods as shown in 
Fig. 3 are not greatly different. 
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APPENDIX (1) 


Including the effect of the heat conduction along the 
skin, the heat flow into a unit surface area of the skin 
(again neglecting radiation) is given by 


Ol ra) oT, 
¢ = G = ; k,,0 : 
Ot Ox Ox 


where k,, is the thermal conductivity of the skin mate 
The more correct version of Eq. (5) is, then (for 


rial. 

kn = constant), 

oT; k, O71, A(t) { d{T,(é, t) — T.(t)] 
- r oo re" x Jt =0 l1 — (¢/x) “*] 


The problem is now vastly more difficult than Eq. (5), 
since, for the step function in velocity, the temperature 
is no longer a function of t/x Some numerical 
integrations of the above equation indicate that for 
most metals the effect of the conduction of heat along 
the skin is very pronounced near the leading edge but 
dies out rapidly away from the leading edge; this is 
to be expected, since large values of 0°7),/Ox? exist 
near the leading edge. 
APPENDIX (2) 
In Eqs. (6) make the following change in variables: 
0 = At/Gs 
Then 
GOr GdTf A 
A oO A dé Gx 
*) (dT/d0) 


= (I/x 


Hence, Eqs. (6) become 


| — (0/e)*7]"" 


dT ‘ P |\d7(&) (dE) dé 
do / 


Let w = d7/dé. .Then, 


w(t) dé 
w(é) = - 
Jo 1 — (0/&) “*| 


Let 0° = 6’ and w(@) = w’(0’). Then, 
— 2 £2 w(t) dé 
WwW v 4 
oO Sa’ \g — @') 
or 
w’ + (2/3)r(2'3)l *w’ = 0 
where 
a I : y 7 + a 
I*f(x) = -— = f(g) (—& — x)*—' dé 
(a Jk 


I*f is closely related to the Riemann-Liouville integral, 
which is the generalization of integration and differen 
tiation to noninteger values. 

Thus, the preceding equation can be interpreted as a 
linear differential equation with constant coefficients 
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Ratio of heat-transfer coefficient to isothermal surface 
coefficient vs. time. Flat plate—exponential 


Fic. 4. 
heat-transfer 
decrease in velocity. 


but with a noninteger derivative (integral). We there- 
fore try the standard technique of substituting an 
exponential function for the dependent variable. 


Then, 
r'(e"™) = —e~™ /*" 


Hence, the characteristic equation of the noninteger 
differential equation is 
1 — (2/3)P(2/3)A~** = 0 
or 
\ = (I'(5/3)] 


Thus, the solution is 


wd) = Ae !P/Sal 


T(0) = A y| e PG/)aP/A dy + B 
0 


and from the boundary conditions B = 7;, A 


T —T, "6 _ 
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Thermal Stress in Power-Producing 
Elements 


A. S. THOMPSON? 
Oak Ridge National Laboratory 


ABSTRACT 


This paper considers the stresses in a material resulting from 
the flow of power through it, when the density of power generated 
has a given distribution. From Hooke’s law relating elastic 
stresses and strains, a general expression is derived for the order 
of magnitude of the thermal stress produced by a temperature 
difference in the material. The special case is studied of the stress 
in a heated cylindrical tube with a given distribution of power 
density and with heat flow through its inner surface. The stress 
in this cylindrical tube resulting from a power distribution in the 
form of a step function of the radius is determined 

It is shown that the expression for the maximum stress in the 
heated tube depends on three groupings of parameters—the 
first involving the properties of the heated material; the second, 
the average power density generated; and the third, geometry 
From these groupings of parameters, the follow- 
(1) A material should be used having 


considerations 
ing conclusions are drawn: 
low values of the modulus of elasticity and the coefficient of ther- 
mal expansion and a high value of thermal conductivity; and 
(2) from geometry considerations the tubes should be of small 
diameter and the power should be generated in a thin layer of 
material directly at the cooled surface. 

These thermal stresses are of importance, since in extreme cases 
they are the cause of spalling (surface failure). 


NOMENCLATURE 


= strain in material 


€ 
ao = stress in material 
E = modulus of elasticity 
u = Poisson’s ratio 
a = coefficient of thermal expansion 
T = temperature in material 
r = radial distance in cylindrical tube 
a = radius of inside of tube 
6 = radius of outside of tube 
q = power density per unit volume of solid material in heater 
6 = average power density per unit volume of the heater 
= thermal conductivity of material 
eo = r/b 
w = a/b 
7, = (T T.)/(Te — Tr) 
y = (q/k)/(q/k-) 
K = Ea(T, T,)/(1 — p) 
Subscripts 
t = tangential direction 
r = radial direction 
z = axial direction 
c = reference point 
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INTRODUCTION 


a BRITTLE MATERIALS exposed to large differences oj 
temperature throughout their volume and in ductik 
materials in which there are repeated cycles of tempera 
ture difference, or surface failure, due to 
thermal stress becomes a significant design problem. 


“spalling,” 


Two classes of problems arise involving thermal stress 
in a material in which power is developed. The first 
consists of those cases where good design can eliminate 
the temperature difference and, hence, the problem 
This category includes temperature variations due to 
nonuniformities in coolant flow, in power density 
throughout the material, and in mechanical structure. 
The second type of problem is inherent in the process 
of transferring heat to a cooling medium. This latter 


problem will be considered in this report. 


ORDER OF MAGNITUDE OF STRESS 


An estimate of the maximum possible elastic stress 
that can exist in a homogeneous material due to a given 
temperature distribution will be made from a considera 
tion of Hooke’s law, by which the elastic strains in three 
mutually perpendicular directions are related to the 


stresses as follows:! 


e = (1/E)[o, — pwloo + o3) + Ea AT 
e = (1/E)loo — pwlos + o;) + EaAT 2 
6, = (1/E)[o3 — wlor + os) + Ea AT 3 


In general, the strains will be different in the three 
directions, and, if strains are due only to temperature 
distribution, they can be expressed by 


¢ = maAdT } 
e = maAdT d 
Eg = nsaAdT 0) 
where 0 < m1, m2, 93, < 1. Here, 7 = 1 corresponds to 


unrestrained expansion due to temperature change, 
AT, and n = 0 corresponds to a complete restriction 0! 
temperature expansion. 


With this expression for strains, the stresses are give! 


by 
a; = -nmkaAdT (@ 
o. = —nmEaAT S) 
o3 = —n3EaAT (9) 
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STRESS IN POWER 





Fic. 1 
where 
| x 
n 
; 1 — 2 
) | { 
| — (1 — wim + w(me + 73) | (10) 
| 1+ pu f 
| 
n x 
| — Pu 
J { 
(1 — w)no + w(x + m1) (11) 
' l+u ( 
l 
nN x 
l Pu 
| { 
) . (1 ns + wm + 12) (12 
| l+-p f 
It can be seen that when yn; = ns = 0, m = me = 
n; = 1/(1 — 2u); and when m = ne = 93 = 1, m1 = M2 = 
n, = 0. Since Eqs. (10), (11), and (12) are monotonic 


in m, 7», and n3, there are no maxima or minima between 
these end values and 


0 < m1, M2, m3 < 1/(1 — 2p) (13) 


The maximum value of the thermal stress which can 
possibly exist, corresponding to full restraint of all 


thermal strains, is, hence, 
2) (14a) 


Can = EeAl./{t = 


is the maximum temperature difference 
The maximum 


where A7’,,,,, 
existing in the homogeneous material. 
thermal stress in a structure will, in general, be less than 
that given by Eq. (14a), its value depending on the con- 
straints actually imposed. 

When there are no restraints in the third direction 
then m; = 0. The maximum stress then occurs when 


m = m = Oand is given by 
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Smear. = EaAT wez./(1 — ps) (14b) 

This is the so-called two-dimensional case of stress. 
Likewise for the one-dimensional case of stress 

Omar. = EadT maz ( l te 


Eqs. (l4a)-(l4e) give the order of magnitude of the 
thermal stress due to a given temperature difference 
in the cases of three-, two-, and one-dimensional stress, 


respectively. 


CALCULATIONS FOR A TUBULAR HEATER 


For any particular case in order to find the actual 
stress at a given point, it is necessary to solve the differ- 
ential equations of the theory of elasticity, for equilib 
rium of forces and compatibility of strains, subject to 
Hooke’s law for the linear range and to the pertinent 
boundary conditions. As an example of this procedure, 
a solution will be made for the stresses due to heat dissi- 


pation into a circular duct. 


Temperature Distribution 


In the core of the heater, heat is generated in the 
material that surrounds the coolant ducts and is dissi- 
pated into the cooling medium flowing through the 
ducts. The idealized situation is shown in Fig. 1. 
Heat is generated in the cross-hatched section at a rate 
g, depending on radial position, and is dissipated on the 
boundary r = a. The boundary 7 = 6 is assumed in- 
sulated. For a steady-state heat balance over a small 
section at r, assuming circumferential symmetry and 
axial uniformity of all variables, it is required that the 
temperature distribution satisfy Poisson’s equation for 
heat flow, which is in cylindrical coordinates for this 


ld (“ ) e q 
r dr ‘@ - k 


Eq. (15a) expressed for convenience in terms of dimen- 


special case.” 


(l5a) 


sionless variables as defined below, is 


d dr b*(q/k), : 
dp (» a . i — 
where 
q = heat added per unit time per unit volume at 
radius, r 

k = thermal conductivity 

7 = temperature at radius, 7 

r = (T — %)/(Ta — To) 

op = r/b 

w = a/b 

y = (q/k)/(q/k)c where c is some reference point 


Integrating Eq. (15b) subject to the conditions 


(16) 


dr dp|, 1 = O 
(17 


T ,= 1 


p w 
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gives, for the temperature distribution in the material, 


p dp, pl 
Y2p2 dp» 
1 Pi 1 


T= . = (1S) 
w dp, pl 
Y2p2 dp» 
1 Pi 1 
where 

“wo dp “pl 7, - Ty 
Y2pP2 dp» = = ( 19) 

l pi l b°(q k). 


Power Density 


The average power density developed per unit volume 
of heater is 


7 


6 = —2q, | G/Qee1 Api 
1 


For the case where the conductivity, k, of the material 
is constant, 


* 


6 = -4. J Yip1 dp, (20) 
1 


Combining Eqs. (19) and (20) to eliminate g, gives, for 
the average power density in the heater tube, 


2k(T, — Ts) | yee 
é=- - ; 
b? * dp, ya 
Y2P2 dp» 
1 pi JI 


Stress in Heated Tube 


(21) 


Because of the temperature variation given in Eq. 
(18), stresses exist in the solid material of which the 
heater is composed. If it is assumed that the cylinder 
is long compared to its diameter, axial strain at points 
away from the ends of the cylinder can be considered 
plane and therefore independent of radial position on 
the cylinder. The equilibrium of forces in radial direc- 
tions for such a cylinder requires that! 

\d(pa,)/dp| — o, = O (22) 
The radial and tangential strains for this case are re 
lated to the deformations by the expressions, where 1 
is the dimensionless deformation. 
e, = du/dp (23) 
€, = U/p (24) 
Eliminating u from these equations gives the com- 
patibility condition for radial and tangential strains. 
»)7 


rAS | 


e, — (d/dp)(pe,) = O 
From Hooke’s law for the elastic case, 
e- = (1/E)[e, — ulo, + a2) + EKa(T, — 7,)r (26) 
e, = (1/E)[o, — w(o, + o,) + Ea(T, — T))r (27) 


e. = (1/E)[o. — ulo, + o,) + Eal(T, — T;)r 


constant (2S 
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Solving Eq. (28) for 
o: = p(o, + o,) + Ee. — EalT, T,)r (29 


Eliminating o, from Eqs. (26) and (29), 


e = (1/E)[(1 — wo, — wl + po, + 
(1 + pw)Ea(T, — T,)r — wke (30) 
e = (I/E)[( — w*)o, — wl + wo, + 


(1+ p)Eal(T, — T;) - uke. 3] 


Substitution of these quantities into Eq. (25) gives the 


compatibility equation in terms of stresses. Since , 


is constant, 
(1 — p)o, — wo, + Eal(T,, — 7),)7r - 
(d dp)}p[(1 — p)o, — wo, + Eal(T, — 74)7]} 0 (3 


Substituting for o, its value from Eq. (22), 


1( Py | d 
(l — p)o;, — “ ‘a + Ka(T, — 7,)r — / 4 
dp dp 
(po, 
a (1 — p) sd — Ur Kall, — 71\,)r Q (33 
| dp f 


Differentiating the products and collecting terms, 


do, , d*a, K dr 
3p + p ~ = ie = 0 34 
dp dp’ dp 


where 
K = Fka(l), — 73)/(1 — pw) 35 


This can be written 


d ( P S:) L Ko? dr 9 oe 
ad = Bin) 
dp r dp - dp 


Integrating Eq. (36) from | to p, 


do l (=) AK ¥ Ld 0 9 
de 2 dp : = 23 JS, Pi i 0 


Integrating the right side by parts, 


do, | /do, K 2K p . 
a ee si a e Tip) dpy =U ” 
dp  p*\dp/y pp p> Si 
Integrating a second time, 
| /do, l , 
CG; (o-)5 + _— | + 
2 \dp/» \p° 
: "p T} : ' dp, “p = 
kK dp, — 2A . Top» dp» = |) aU 
Ji pl 1 pi” JI 


. 


! 


+t 


Integrating the last member by parts and collecting 


(=) ( ) 
i fas 
2 \dp/» \p 
A p 
Tip, dp, Q (40 
p” Ji 


Since no radial constraints act at the boundary, 


terms, 


G; (o,), + 
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(g-). = (es = 0 (41) 


Substituting into Eq. (40), 


es K f° 
(“) ( ') + : / 71p1 dp, = 0 (42) 
) dp » \@ wo Ji 


Hence, 
14 kK : 
? ( dp h 3 l : w" 


/ Tipt d py (45) 
é J | 


The expression for radial stress now becomes 


. 


K | p- = 
0, = ( . / TiP\ dp, 
p \ Ww | af 1 
/ Tip dp.) (44 
rt 


From Eq. (45), 


do 2K o K 
— “7s Tip; dp, — ie 
ga" — | J 1 p 


> 
p 


2h Of 
J Tip: dp, (45 
pr Ji 

Substituting these values in Eq. (22) for tangential 
stress, 


K f 
a, = ( Tip) dp, T 
Pp \/1 
1 + p- *w 
= , Tip! dp, — Tp (46 
ee 1 


w . 


From Eq. (29) the axial stress is 


: 2u ’ 
a K ( / T1p; dp, — r) + le (47 
Ww wt 


Since axial strain is plane and there are no end con 
straints, 


[o-1do = 0 (4S 
rl 
2 “a 
¢; = K ( / Tip: dp, — r) (49 
Ww Pe 


By equating da,/dp 


Hence, 


radial stress is shown to be 


Os wie (A 2)| fw Fd (50 
where 
. ») ws 
I(w = / TiP1 dp) ol 
w L ova 
Hence, the total possible variation of o, mar, 1S 
(K/2)| fw li S thas. S (2/2 (52 
rhe greatest tensile stress exists at 7 a. Itis 
ae 0 K | f(w) | (53 
rhe greatest COMpressi\ e stress exists at ? b. Itis 
Co; 0 Af(w (o4 


0 the maximum value of the 


It can be seen that the maximum tangential and axial 
stresses will always be of the order of twice the maxi- 
mum radial stresses. These stresses are in each case 
the product of K, the order of magnitude term as given 
by Eq. (15b) for the case of a two-dimensional stress 
problem and a quantity that is a function only of the 
geometry of the particular tube under consideration. 
To evaluate these stresses, it is necessary to evaluate 


{(w). Substituting the value of 7 from Eq. (18) 


i ‘do [* 
pr dp) 733 dps , 
Ji 1 Pp2 J1 3 
“ede 
| / Y2p2 dp» 
l Pi J! 


(a) = 


(55 
Integrating the numerator once by parts gives 
| pi dp, / 72p2 Up» ; 
1 Jil 
Mw nas w” = . . > 
” dp, 4 Ww _ l 
y2p2 dpe 
l Pi Ji 
(56) 


Integrating both numerator and denominator by parts, 


ri wo? — 
J\ 2 - 


Yip! dp) 


Ww 7 aw Foy - gt l 
log — yipi dp; 
J 1 Pi 


Assume a radial distribution of power generation in the 
material around a cooling hole as shown in Fig. 2. It 
can be seen from the figure that heat is generated at a 


uniform rate per unit volume in the region w < p < 8B. 
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No heat is generated in the region 8 < p < 1, which is 


therefore maintained at a constant temperature 
throughout. For this case, 
w w" ares, a 
Pi dp, 
: 8 2 7 
fw) = | w? — - - (58) 
ss @ = 
| log — p; dp; 
B Pi 
Performing the indicated integrations, 
wi — Bt — 4w*B* log (w/ B) , 
tlw) = (59) 


2(w? — 1)[w? — B? 


2B" log (w/ B) 


By taking limits as w varies between 0 and 1, it is found 
that 


Os f(w) Ss! 


Oo 


IIA 


] (60) 


oy 
or 


*/, S 1 — f(w) = J (61) 


Hence, it is seen from Eqs. (53) and (54) that the maxi- 
mum stress in the heater is tensile and has a value, at 
r = a, between the limits 


IIA 


K (62 


6 3K > (Ot mar.» Cz maz.) ) 


regardless of the configuration of the heater tube and 
regardless of the extent of the heated area. The only 
hope for material reduction of temperature stress is 


then to reduce A. From Eqs. (21) and (35), eliminating 


(Ts Se Fs); 
K Kab*6 | 2 | w 6 
= — - og (63) 
(1 — wk (w/s)?-1 "pl *‘ 
This expression goes to zero for w/8 = 1. Hence, fora 


given power density, it is desirable to concentrate the 
heating effect as close as possible to the cooled surface 
of the heater tube, this giving the lowest possible stress 
due to dissipation of a given power density in the tube. 
From Eq. (63) it can be seen that, for a given tube con- 
figuration and depth of heating, the maximum stress 
imposed on the heater by withdrawal of power is regu- 
lated by the magnitude of the factor Eab?6/(1 — w)k. 
Hence, for a given power density it is desirable to have 
small heating tubes in a material having low values of 
modulus of elasticity and coefficient of thermal expan- 
sion and a high value of thermal conductivity. 

It is to be noted that the order of magnitude of the 
thermal stress is found in Eq. (15a) to depend primarily 
on a temperature difference (the integral of temperature 
gradient over a path) rather than primarily on tempera- 
ture gradient as maintained by some authors.* Thus, 
the stress due to an infinite temperature gradient, if the 
temperature difference is finite, is not infinite but is 
merely indeterminate from the point of view of temper- 
ature gradient. 


Example 


Consider for an example a heater having the following 
characteristics: 





5b = 0.5.cm. 

w = 0.554 

B = 0.612 

6 = | kw. per ce. 

Ek = 20,000,000 Ibs. per sq. in. 
k = 0.5 watt cm./cm.? °C. 

nh = 0.5 

a = 20 X 10° per ~C. 


Substituting in Eq. (59), the value of /(w) is 0.03. The 
quantity [1 — f(w)] is, hence, 0.97, or for all practical 
purposes 1, so that the maximum stress is given by the 
value of —A. Substituting values in Eq. (63) gives 
Omar. = —K = +6,400 Ibs. per sq.in. (tension), 


SUMMARY 


Eqs. (l4a)—(l4e) are general equations for the order 
of magnitude of thermal stress associated with heat 
transfer. For the special case of a heated cylindrical 
tube, with a given distribution of power density and 
with heat flow through its inner surface, Eqs. (50), 
(53), and (54) show the maximum (tensile) and mini- 
mum (compressive) radial, tangential, and axial stresses 
as the product of an order of magnitude term, A, anda 
geometry factor. The maximum radial stress is of the 
order of one-half the maximum tangential and axial 
stresses. The latter are, for all practical purposes, 
equal to the order of magnitude term, A, since the 
geometry factor is nearly unity. 

‘The evaluation of the maximum thermal stress hence 
consists in evaluating A. Eq. (63) shows this to de 
pend on three groupings of parameters —the first in 
volving the properties of the heated material Ka 
(1 — w)k; the second the average power density 6; and 
the third geometry considerations 


2 w 
b? (: — log ) 
(w/ Bp) — 1 B 


The factor ) in the last expression is a measure of the 
size of the tube, and the parenthesis involves the di 
Neither 
the temperature difference nor the temperature gradi 


mensionless shape of the power density curve. 


ent appears explicitly in Eq. (63), although the tem 
perature difference is implicit in power density consider- 
ations. 

To minimize the stresses in a material due to a given 
power density, it is therefore desirable to choose a ma 
terial with a low value of Ea/(1 — w)k, to use small 
tubes (small 5), and then to develop the power in a 
thin layer of material directly at the surface (w = 38). 
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Visualization of Flow Fields by Use of a Tuft 
Grid Technique 


JOHN D. BIRD* 
Langley Aeronautical Laboratory, N.A.C.A. 


SUMMARY 


A technique for obtaining ‘a physical picture of the flow field 
behind a wing, combination of wings, or other airplane com- 
ponents is described wherein the action of a great number of 
tufts of uniform length, attached to a wire grid, are photographed 
from far downstream. This procedure yields, with a minimum 
of labor, an approximate vector plot of the flow field in a plane 
normal to the air stream. An examination of the data obtained 
from an application of this tuft grid technique to the study of the 
flow fields behind several low aspect ratio rectangular and delta 
wing arrangements indicates that analyses of such factors as 
trailing vortex strength and location and approximate downwash 
and sidewash angles over a large area around lifting surfaces and 
other aerodynamic forms may be satisfactorily conducted by this 


technique 
INTRODUCTION 


— TRENDS IN AERODYNAMIC DESIGN are to- 
ward the use of unorthodox wing plan forms, 
ratios, and unusual combinations of 
lifting surfaces. This is particularly true in the case 
of missile design where low aspect ratio cruciform 


lower aspect 


arrangements and tandem surfaces of approximately 
equal area are commonly encountered. These factors 
require a good understanding of the flow field existing 
in the vicinity, and downstream, of these lifting sur- 
face arrangements in order to provide for an optimum 
choice of geometry and, in particular, for optimum 
placement of rear surfaces relative to those in front. 

A technique for obtaining a physical picture of the 
field kehind a wing, combination of wings, or other 
airplane components has recently been developed at 
the Langley Aeronautical Laboratory, N.A.C.A. This 
technique involves photographing the action of a 
great number of tufts of uniform length, mounted on a 
screen, from far downstream. This permits, in an 
approximate fashion, obtaining a major end result of a 
pitot-yawhead survey—namely, a vector plot of the 
flow field in a plane normal to the air stream at a station 
in the wake of an aerodynamic surface. In this sense, 
a tuft grid should be particularly useful in instruction 
in aerodynamics 

Tufts have, of course, frequently been employed as 
an aid to visualizing flow in the vicinity of aero- 
dynamic shapes and in various forms of wind channel. 
Papers by Purser, Spearman, and Bates! and by Schlicht 
ing* describe the use of tufts in recording the downwash 

Presented at the Aerodynamics Session, Twentieth Annual 
Meeting, 1.A.S., New York, January 28-31, February 1, 1952 

* Head, Stability Tunnel Section 


and general flow character behind lifting surfaces by 
making photographs of tufts suspended on _ wires. 
These photographs were made along a line about normal 
to the free stream and are an excellent means for 
establishing approximate downwash angles with a 
minimum of labor. 


DESCRIPTION OF TUFT GRID AND TEST TECHNIQUE 


The tuft grid used in the application discussed herein 
consisted essentially of a rectangular grid of fine wires 
of 1-in mesh supported at the periphery by a tubular 
framework with 3-in. woolen tufts attached at 1-in. 
intervals over the surface of the grid. The tufted area 
of the grid covered an area 50 in. wide and 26 in. high. 
A preloading system consisting of a spring mounted 
between one end of each wire and the frame was used 
to maintain a 1.5-lb. tension in each wire. Each inter- 
section of a vertical and a horizontal wire was soldered 
so that relative movement between wires would be 
eliminated. 

Three-inch lengths of four-ply wool knitting yarn 
were used as tufts, with single tufts attached to each 
intersection of a vertical and horizontal wire with strong 
thin twine. A small loop of twine was provided at the 
attaching point in order to permit the tuft to move 
freely in all directions, and the downstream end of each 
tuft was tied with thread to prevent the strands of 


wool from unraveling. 


A trusslike structure was welded to each of the 
rectangular frame members to provide the necessary 
rigidity. The frame and truss structure were fashioned 
from thin-walled streamline tubing having a 0.75-in. 
chord. Mounting brackets were located at each of the 
four corners for attaching the frame to the supporting 
system used in the tunnel. The supporting system 
was designed so that the screen could be moved up- 


or downstream at will while the tunnel was in operation 


A cutaway view of a portion of a wind tunnel showing 
the general location of the camera, grid, and model is 
presented in Fig. 1. The camera, of course, can onlv 
record the projections of the tufts and model in a ver- 
tical plane. Fig. 2 gives an illustrative example of 
the type of picture obtained by a still camera. In 
this figure, the gray triangular area represents the 
vertical projection of a triangular-wing model at an 


angle of attack. The short dark lines represent the 
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lic. 1. General arrangement of model, tuft grid, and camera for 
testing. 
. 
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Fic. 2. Typical flow pattern obtained by tuft grid technique 
Triangular wing 
tufts. In the relatively undisturbed regions of flow, 


the tufts appear as short dark lines or dots, as may be 
seen at the center of the top horizontal row and in 
the disturbed regions as longer lines. One of the tufts 
in the sketch has been magnified for the purpose of 
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illustration. The vertical and horizontal Projection 
of the tuft is an indication of the local angles of down. 
wash and sidewash, respectively. The vertical] and 
horizontal projections of each tuft, together with the 
tuft length, which is 3 in. for this grid, determine the 
downwash and sidewash angles of the flow at a given 
location. A parallax correction proportional to the 
distance of an individual tuft from the center of the 
screen must be applied to measured angles. Deter. 
mination of the change in downwash at a specific 
point caused by some minor model modification sych 
as filleting could better be accomplished, of course, by 
the use of a yawhead. A practical limit of accuracy 
of the tuft grid employed in this investigation for 
regions of reasonably steady flow is felt to be +! 
of angle for an individual tuft. 

Some efforts have been made by Mr. Boatright, oj 
the Langley Aeronautical Laboratory, N.A.C.A. to 
develop a tuft grid suitable for use at supersonic speeds 
His chief difficulty to date has been in obtaining a tuft 
that is sufficiently stable for satisfactory testing at 


these speeds. 


RESULTS AND DISCUSSION 
Presentation 


An investigation of the tuft grid technique as de 
scribed herein has recently been completed, and some 
of the results of the flow field studies behind a number 
of arrangements of aerodynamic surfaces are presented 
to demonstrate the range of applicability of this tech 
nique. Tuft grid surveys pertaining to the flow fields 
behind rectangular and 60° delta wings, the generation 
of the trailing vortex system across the chord of a 
rectangular wing of '/; aspect ratio, and the occurrence 
of the leap-frogging phenomena behind a lifting X 
wing configuration are used as examples. 


Surveys Behind Rectangular and 60° Delta Wings 


Results of a tuft grid study of the flow behind a 
rectangular wing of 2.61 aspect ratio and a 60° delta 
wing of 2.31 aspect ratio are shown in Figs. 3 and 4 
The tuft grid was located 
2 ft. from the model in each case. The rectangular 
wing was of 3-ft. span and had an NACA 0012 airfoil 


for several angles of attack. 


section. The delta wing was also of 3-ft span and hada 
6 per cent thick airfoil section in a streamwise direction. 
The effect of sideslip on the flow field of the 60° delta 
wing is shown in Fig. 5. 

The increase in strength of the trailing vortex systems 
of the wings with increase in angle of attack may be 
noted in Figs. 3 and 4. The trailing vortices for the 
rectangular wing remain essentially at the wing tp 
throughout the angle-of-attack range, and the stall 
becomes evident at the center of the wing at 1S 

The trailing vortex system for the triangular wings 
moves inboard as the angle of attack is increased, how- 
ever, indicating the formation of vortex trails above 
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VISUALIZATION 


the surface of the wing as discussed in a paper by 
Wick,’ of the Ames Aeronautical Laboratory, N.A.C.A. 
These vortex trails are a result of a chordwise rolling-up 
of the trailing vorticity and may occur for swept or 
delta wings. The trails develop at the apex of the 
wing in conjunction with a slight leading-edge separa- 
tion and continue across the chord slightly above the 
surface of the wing. The regions of the wing outboard 
of these trails are stalled, of course, and carry little 
load. The stalled region increases in magnitude and 
the trails move inboard as the angle of attack is in- 
creased. 

A comparison of the spanwise positions of the trailing 
vortices of the rectangular and delta plan-form wings 
studied for this investigation are given in Fig. 6 as a 
function of angle of attack. These data were obtained 
from the tuft grid photographs of Figs. 3 and 4. The 
inward movement of the trailing vortex system of the 
delta wing is seen to be appreciably larger than that of 
the rectangular wing. The effects of sideslip on the 
flow’field behind the delta wing are to make the trailing 
vortex from the retreating panel more steady and more 
tightly rolled up and to cause the trailing vortex from 
the advancing wing to become unsteady and to exhibit 
the characteristics of an unswept wing approaching 
the stall (see Fig. 5). 

These surveys emphasize the fact that horizontal and 
vertical tail location will have more influence on the 
contribution of tail surfaces to longitudinal and direc- 
tional stability for the low aspect ratio wings that are 
increasingly gaining favor than has been the case for 
higher aspect ratio wings. Variations in angle of 
attack and sideslip will bring tail surfaces into, and out 
of, the powerful effect of the trailing vortex system in 
direct dependence on the placement of the tail with 
respect to the wing-chord plane 


Surveys Along Chord of Rectangular Wing 


The results of an application of the tuft grid technique 
to the study of the development of the trailing vortex 
system along the chord of a rectangular wing of !/; 


aspect ratio having a flat-plate airfoil section are given 





Fic. 3. Tuft grid surveys 2 ft. behind a rectangular wing of 
aspect ratio 2.61 having an NACA 0012 airfoil section for angles 
of attack of 0°, 12°, 16°, and 18°. Wing span, 3 ft.; dynamic 
pressure, 8 Ibs. per sq.ft. 
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Fic. 4. Tuft grid surveys 2 ft. behind a 60° delta wing having 
a 6 per cent thick airfoil section for angles of attack of 0°, 12°, 
16°, and 20°. Wing span, 3 ft.; dynamic pressure, 8 Ibs. per 
sq.ft. 





delta wing at 20 
dynamic 


Tuft grid survey 2 ft. behind a 60 


Fic. 5. 
angle of attack and —10 
pressure, 8 lbs. per sq.ft. 


sideslip. Wing span, 3 ft 











WING TIP ae 
TRAILING RECTANGULAR WING 
Boation 60° DELTA WING - 
WING ROOT t : 
0 10 20 
a, DEG 
Fic. 6. Spanwise locations of trailing vortices 2 ft. behind 


rectangular and 60° delta plan-form wings of 2.61 and 2.31 aspect 
ratio, respectively, as a function of angle of attack. Data ex 
tracted from Figs. 3 and 4. 


in Figs. 7 and 8. For this study, a slot equal to the 
span of the wing was cut in a tuft grid to permit mount- 
ing the grid at various chordwise positions on the wing. 
Tuft grid surveys at chordwise stations of 0, 12'/s, 
2, 00, and 75 per cent of the wing chord are 
presented for a 20° angle of attack at a dynamic pres- 
The progressive rolling- 


25, 37! 
sure of 8 Ibs. per sq. ft. in Fig. 7. 
up, increase in strength, and vertical displacement of 
the trailing vortex system relative to the wing-chord 
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Fic. 7. Tuft grid surveys at chordwise stations of 0, 12'/2, 25, 
37'/2, 50, and 75 per cent air foil chord on a rectangular wing of 
'/, aspect ratio. Flat-plate airfoil section; dynamic pressure, 8 
Ibs. per sq.ft. 





VALUE CORRESPONDING TO 
MEASURED LIFT—— 








Teese. ff  — | SEEReaeanse en 
VORTEX 
STRENGTH 
t 
LEADING TRAILING 
EDGE EDGE 
Fic. 8. Strengti of trailing vortices as a function of chordwise 


station for a rectangular wing of '/; aspect ratio. Data obtained 


by contour integrations from Fig. 7 


plane are clearly shown in these surveys. The down 
ward slope of the trailing vortex cores determined from 
these photographs is 12'/,°. This value is the same as 
the value predicted by the theory of Bollay.' 

The results of contour integrations of the tangential 
velocity around the individual trailing vortices showing 
the increase in trailing vortex strength across the chord 
of the wing are presented in Fig. 8. These contour 
integrations were made from measurements of the 
projections of tufts from the photographs in Fig. 7. 
The projections were used to determine the required 
sidewash or downwash angles (see Fig. 2). The agree 
ment between the trailing vortex strength obtained at 
the trailing edge by contour integration and the trailing 
vortex strength corresponding to the measured liit is 
surprisingly good. 

The curve of vortex strength is, of course, approxi 
mately proportional to the accumulation of load across 
the chord of the wing. In this sense, the gradual 


build-up of load across the chord as indicated by these 
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results is markedly different from the linear theory 
of Wieghardt,® which predicts a much greater cop. 


centration of load near the leading edge. The center 


of pressure as predicted by linear theory is 7 per cent 
of the chord from the leading edge as contrasted to the 
results of force tests on this '/4 aspect ratio wing which 
indicate the center of pressure to be at 36 per cent 
of the chord for a 20° angle of attack. The lift pre. 
dicted by the linear theory is, of course, only about 


one-third of the experimental value. 


Surveys Behind an X Wing 


A phenomenon that is closely related to a similar 
behavior of vortex rings as described by Helmholtz 
in 1858 is illustrated by the tut grid surveys shown in 
Fig. 9. This phenomenon has been called leap-frog- 
ging and is associated with the behavior of the trailing 
vortex system of an X arrangement of lifting surfaces. 
The leap-frog distance is defined as the distance from 
the wing trailing edge to the point downstream at 
which the trailing vortices from one pair of tips Have 
passed inside the trailing vortices from the other pair 
of tips. This distance is dependent on the aspect 


ratio and lift of the surfaces Tuft grid surveys at 





Fic. 9 2-, and 4-ft. stations behin¢ 


Tuft grid surveys at ! ; 
an X arrangement of 70° triangular wings at 20° angle of attack 


Flat-plate airfoil section; wing span, 0.8 ft.; dynamic pressure, 


8 Ibs. per sq.ft. 
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VISUALIZATION 


|, 2 and 4 ft. behind the trailing edge of an X arrange- 


ment of triangular surfaces with 70° of sweep of the 
leading edges are shown for an angle of attack of 20° 
and a dynamic pressure of § Ibs. per sq.ft. in Fig. 9. 
The wing span was 0.8 ft. White strips of tape were 
attached to make the trailing edges of the wing 
visible. 

Four discrete vortex cores may be seen at the !/2- 
ft. station. The upper pair is considerably inboard of 
the lower pair 
station where the upper pair of vortices is no longer 
easily discernible because of proximity to the lower 
pair. At the 4-it. station, leap-frogging is complete, 
and the upper pair of vortices has passed through and 
The appearance of the flow 


Leap-frogging occurs near the 2-ft. 


beyond the lower pair. 
field at this station indicates that some dissipation or 
spreading of the inner pair of vortices has taken place 
during leap-irogging. 

The general appearance of these surveys indicates 
that the placement of tail surfaces behind a low aspect 
ratio X wing would require special consideration of the 
trailing vortex system. The strong downwash fields 
that shift nonlinearly with angle of attack and sideslip 
may introduce undesirable variations in the contri- 
bution of the tail surfaces to the aerodynamic forces 
and moments necessitating unusual variations of con- 
trol surface settings for trim. 
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CONCLUDING REMARKS 


An examination of the data obtained from an appli- 
cation of the tuft grid technique to the study of the 
fluw fields behind several low aspect ratio rectangular 
and delta wing arrangements indicates that analyses 
of such factors as trailing vortex strength and location 
and approximate downwash and sidewash angies over a 
large area around lifting surfaces and other aerodynamic 
forms may be satisfactorily conducted by this tech- 
nique. In addition, excellent illustrations of aero- 
dynamic phenomena involving the interaction of the 
combinations of lifting surfaces are 


flow fields of 


possible 
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Buckling of Sandwich Cylinders Under Bending and 
Combined Bending and Axial Compression 


(Continued from page 470) 


( a) C. + ae ae + Cru) = 0 
Q n IC n—1 n+l _ 


or, by dividing through the common factor .\/at,7C, we 


find 
Cat L'C, + Can = 0 (17) 
for m * 1, where L’ = [1 — (N,/C)] IC/Mat. For 
n = 0), we have 
L'G + CG, = 0 (18) 
and, form = 1, we have 
20% + L'C, + CG. = 0 (19) 


Eqs. (17), (18), and (19) are exactly the same as Eqs. 
(9), (10), and (11) if we replace L’ by L. The solution 
istherefore L’ = +2 or 


N, + (2Mat/I) = C = hG (20) 


Eq. (20) again indicates that buckling of sandwich 
cylinder will occur when the maximum compressive 


stress in the cylinder is equal to the buckling stress when 
the cylinder is under compression alone. 
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On Pohlhausen’s Method with Application 
to a Swirl Problem of Tavlor 


J. C. COOKE* 
Unwersity of Malaya 


SUMMARY 


Most applications of Pohlhausen’s method to three-dimensional 
laminar boundary-layer flow have assumed a ‘‘thickness of the 
’ 


boundary layer,”’ 5, at which the velocity in the boundary layer 


becomes equal to the velocity in the main stream. However, in 
solving problems, the velocities over the surface are split into 
components at right angles, and each component is supposed to 
reach the corresponding main stream component in the same 
distance 6. It is asserted in this paper that this is not correct, 
and examples are given to show that the assumption leads to re 
sults much less accurate than those of which the method is ca 
The method is illustrated by the solution of a swirl prob- 


The algebraical and numerical work is 


pable 
lem solved by Taylor. 
heavy even for cases where there is symmetry of some kind, as 
there is in Taylor’s problem. 


INTRODUCTION 


| pew EXAMPLES EXIST of the solution of three-dimen- 
boundary-layer Sears,! Wild,’ 
Young and Booth,* and Cooke* * have given solutions 


sional problems. 


for infinite yvawed cylinders; Taylor’ and Binnie and 
Harris’ have given solutions of swirling problems; 
Carrier’ has given the solution of the problem of the 
boundary layer in a corner. There are other solu- 
tions,?~!? one of which" is an exact solution of the full 
Navier-Stokes equations, and the others are exact 
Pohl- 


hausen’s method is the only one to have been used in 


solutions of the boundary-layer equations. 
the approximate solutions, and Cooke® found that this 
method is capable of good results in one series of in- 
stances when the pressure gradient is not adverse. 
All those using the method except Wild? and Cooke? 
assumed that the two components of velocity in the 
boundary layer reached the corresponding component 
in the main stream in the same distance 6, the ‘‘thick- 
ness of the boundary layer.’’ Sometimes this leads to 
incorrect results, though with luck in the selection of 
the polynomials it may be satisfactory. 


THE BASIS OF POHLHAUSEN’S METHOD 


We consider first the steady two-dimensional lami- 
In the 
many exact solutions obtained, it is always found that 


nar flow of a fluid stream close to a fixed body. 


the velocity in the boundary layer approaches asymp- 
totically that of the main stream, l’, and so it is not 
possible to specify the boundary-layer thickness ex- 
actly, but the layer is taken as extending from the body 
Received April 2, 1951. 
* Professor of Applied Mathematics. 
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to the point where the velocity has reached, say, 0.995 
that of the main stream. 

Pohlhausen’s method involves assuming a definit, 
thickness, 6, for the boundary layer; the velocity js 
then supposed to be equal to L’ when at distance § from 
the body. Moreover, the velocity is assumed to be 
expressible as a polynomial in z, the distance from the 
body, and the polynomial is chosen to make the velocity 
“merge’’ into U’ as “asymptotically” as possible at z = 
6. The higher the degree of the polynomial, the better 
this can be done but the more complicated the work be- 
comes. Pohlhausen used a polynomial of the fourth 
degree; others have used higher degrees up to the eley- 
enth. The method can be made to work well, espe- 
cially in the cases when Ll” does not decrease as we go 
downstream. 

When we try to extend this to flow over three- 
dimensional surfaces we run into difficulties. Compu- 
tationally, it all becomes more complicated, but still 
possible. In addition, however, it is difficult to de- 
fine 6, because, although the resultant velocity must 
still approach, say, 0.995 that of the mainstream in 
some definite distance 6, on going into details we must 
split the velocity parallel to the surface into compo 
nents, and these do not reach the required percentage value 
of the main-stream components at the same distance. 
Taylor, Binnie and Harris, and Young and Booth as 
sume that they do, but this is not correct. Consider, 
for instance, the case solved exactly in reference 4 for 
Bg = 1. 
that of the main stream at about Y = 2.6 (Hartree);' 
in the other, at about Y = 3.3 (Cooke),* where ¥ is the 


In one direction the velocity reaches 0.995 


distance from the body measured in the same non- 
dimensional units. The position where the resu/tant 
velocity is 0.995 that of the resultant main stream varies 
according to the relative magnitude of the two main 
stream components, but it lies somewhere between 
Y = 2.6 and Y = 3.3 and is at about Y = 3.1 if these 
components are equal. 

It might be said that Pohlhausen’s method is only 
approximate anyway and that, in fact, 6 turns out to 
be a somewhat vague parameter which varies accord- 
ing to the degree of the polynomial chosen and that no 
serious error will result by assuming the ratio of the two 
Actually, the very fact that 6 does 
Thus, 


6’s to be unity. 
vary in this way must be taken into account. 
in the two-dimensional flow over a flat plate we have 


6(U/vx) ? = 3.46, 5.48, 4.64, or 5.83 according as the 
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polynomial is of the first, second, third, or fourth de- 
orees. It still does not necessarily follow that all we 
oid do is to take polynomials of equal degrees for the 
fow to have equal 6’s; in fact, it seems unlikely. It 
must be admitted, however, that we may be lucky in 
the choice of degrees for the polynomials so that the 
js do become very nearly equal. In this case our solu- 
tion is bound to be more accurate than the one assum- 
ing different 6's. This is because (having one less un- 
known) we can introduce an extra parameter into one 
of the polynomials, thereby refining the method. Un- 
fortunately, there is no means of determining a priori 
or a posteriori what the degree of the polynomials 
should be. It appears then that it may be possible by 
luck to obtain a more accurate answer by assuming 
equal 6's, but it would be safer to assume them differ- 
ent. 

One series of instances may be tested by a known 
This is the case of the yawed infinite 
Here, the 


exact solution. 
cylinder considered in references 4 and 5. 
main-stream velocity normal to the generators is taken 
as U = cx™ and along the generators as a constant IV’. 
Reference 5 assumes different 6’s and obtains good re 
sults for accelerated flow with quartic polynomials. 
If, however, we may assume legitimately that the 6’s 
are equal, we may introduce an extra parameter a into 
the v velocity distribution. (The u velocity distribu- 
tion is already known to be quite good as it is and so 
need not be altered.) If we take 8B = 2m/(m + 1) 
and use polynomials of the fourth degree of the form 


u/U = 2n — 2n* + nf + (1/6)A(1 — n)? 
v/V = 4n* — 3n1 + a(n — 3n* + 2n') 
n = 2/6 


satisfying as many as possible of the usual boundary 
conditions, we obtain results as shown in Table 1, 
where the nondimensional skin friction, displacement 
thickness, and momentum thickness in the v direction 
are shown for three values of 8. The notation is that 
of reference 5. 

We cannot expect a good result for 8 = —0.1988, 
since the polynomial for u/U is known not to be good 
in this case. It 7s, however, known to be good in the 
cases 8 = 1 and B = 2. 

On the other hand, if we try out this method on the 
first case on which it was ever tried (namely, the case of 
rotating disc dealt with by von Karman’ and later by 
Cochran"), we find that the solution with unequal 6's 
is not quite so good as that as given by von Karman® 
as corrected by Cochran") where equal 6’s were used. 
It is significant perhaps, that, as Cochran pointed out, 
only polynomials of certain degrees gave solutions that 
were physically acceptable. It would indeed be a 
happy chance if only the combination of polynomials 
that give a physically acceptable solution is the one 
that is the most accurate, but unfortunate’y it is not 
true, as the examples in Table 1 show. 


TABLE | 
oe |: ig i. 
3 = 1.0 Equal 6’s 0.4049 1.1650 0.3833 
Different 6's 0.5773 1.0393 0.4044 
Exact? 0.5704 1.0265 0.4069 
3 = 2.0 Equal 6’s 0.3410 1.1623 (). 3524 
Different 6's 0.6099 0.9828 0.3872 
Exact 0.6053 0.9759 0.3852 
3 = —().1988 Equal 6’s 0.5126 0.8235 0.0423 
Different 6's 0.3968 1.512 0.5920 
Exact 0.3258 1.6696 0.5937 


We proceed to illustrate the method by solving 


Taylor’s problem.*® 


THE SWIRL PROBLEM OF TAYLOR 


“Liquid enters tangentially a chamber of radius R; 
and emerges through an orifice of radius R. which is 
considerably less than Rs after passing down a con- 
verging cone whose vertex subtends an angle 2a... 
An attempt is made to solve a simplified problem in 
which the longitudinal or axial component of velocity 
outside the boundary layer is neglected, only the swirl 
being considered. The main body of the liquid is, in 
fact, taken to be moving with velocity Q/r in circles 
round the axis, where r is the distance from the axis, 
and Qis the circulation.”’ 

We use spherical polar coordinates with the origin 
at the vertex of the cone, 6 = 0 being its axis, about 
which the motion is symmetric. All derivatives with 
respect to ¢ are therefore zero. 

The general equations of motion are given in Gold 
stein'4 whose notation we use; and when the usual 
approximations for the flow in the boundary layer close 


to the cone are made, we obtain 
lop ov Ou 
po OR R* of 


Ou v Ou w- 
uu - - 
OR R 06 R 


1 1 Op 
0 (2 
p R O60 
Ow v Ow wu vy Ow 
u i 2 = ‘ : 3 
OR R0O# R R? 06 
together with the equation of continuity 
Ou 2u 1 ov 
= U0 | 


oR ' RRO 


Eq. (2) shows that variations in p throughout the 
thickness of the boundary layer may be neglected, so 
that p is the same inside the layer as it is in the main 
stream just outside it. 

Outside the boundary layer, 

w=Q R sin 6 
and so the pressure is given by 


Pp l 0? 
=— —————=-_ + COUSTAII 
p 2 R? sin? 0 








-* 


; ~< - FEaee 
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and Eq. (1) becoines 


Ou v Ou 2? w* vy Ou ; 

“OR RW Risin?@ RR? ~~ 

Taylor shows that, if variations in sin @ in the bound- 

ary layer are neglected and sin @ is taken as sin a, the 

integrals across the boundary layer of Eqs. (la) and 

(3), taking account of Eq. (4) and of the fact that u 

vanishes on the cone and outside the boundary layer, 
while v vanishes on the cone, may be written 


2f +> +2 fap 
a ¢ é ( 
oR De hal 


/ 0° Ww “py O'N 
( ee eae ees ) dé = / 16 
/ R?® sin® a R J OF 


and 


- [ (RoR + 20) do + 
- 21 
R? sin @ ¥ OR A 


© @ (uw) do +3 * Uw sn ~» Ow Po 
— 2s” J» de’ 


where the integrals are taken across the boundary layer 
that is, from 6 = a to some other value to be con- 


sidered later. 


THE POHLHAUSEN POLYNOMIALS 


We shall take as boundary conditions for « and w 
the following: 


O"u . 
— = ()?/yR sin? a 
O06" 
2/62 = 0 at@ =a (5) 


até = a — (6/R) 


at 0 = =f 
eee 


and A being the ‘“‘thickness of the boundary layer’’ 


u=0, Ou/d0d = 0 


sin a 06 


~ 


in the « and w directions, respectively; 6 and A are 
functions of R only. 
We write 


0262 
i= = a—————- (=~ ay > g) = 
tyR® sin? a 
026? 
rere AD, 
vR® sin? a 
(6) 
Q wr 
yi => - —_— = 
Rsin a \2 2 
Q 
ane : 7 c) 
R sin a 


where 


n = R(a — 0)/6 and ¢ = R(a — 8)/A. 


Here, we differ from Taylor in two respects: (1) Of 
course, 6 and A are different; and (2) since we have 
introduced an extra unknown into the equations, to 
counteract this we make use of an extra boundary 
condition in Eq. (5)—-namely, that 


07u/06? = Q?/rR sin? a 


at the wall. This condition arises from Eq. (1a), sine 
u =v = w = Oat the wall. Taylor's coefficient of 
{(n) was written Q//R sin a, where was an unknowp 
function of R to be determined. On trying out Tay. 
lor’s solution for FE, it will be seen that this boundary 
condition is not fulfilled. Of course, the polynomials 
never do fulfill all the boundary conditions. We haye 
also raised the degree of g to make it of the same 
degree asf. This enables us to fulfill one more bound- 
ary condition than Taylor, viz., 0°w/00? = 0 at the 
wall which follows from Eq. (3). 

We substitute Eqs. (6) into Eqs. (la) and (3), and 
we also write 

28° R 


A = ; : ; r= 
Ro*?v sin a Ry 


A = Ké, n= KE 
where RX, is the value of R at the circle of entry of the 
fluid. A will be a function of 7 only. 

. After some straightforward but intricate algebra, we 
obtain 


Ar’ d2 of 
(2a — b) + (26 — 4a) +c = - 7 
a ‘ On/ , =0 
AN’ Le 
(net 5dte- 5h 5! + 
d° K’ yl 1 /dg 
ita SOS SS ew x i 
ia K r* rk \o¢ P 


dashes denoting differentiation with respect to r. 
In these equations 


a = IF dy, b= S nf’ dy, c= rg (1 — g’) dn 
d= f nf'dn, e= Sf fdn, k= S fe dn 
j= f cfe’ dn, h = S nf'g dn 


dashes here denoting differentiation with respect to 
the argument of the particular function /(m) or g(f 
concerned. 


THE LIMITS OF THE INTEGRALS 


In the first place, the Eqs. (6) hold only for 7 or ¢ less 
than unity. As soon as 7 reaches unity, u becomes 
zero and remains at zero for greater values; while, 
when ¢ reaches unity, g(¢) becomes equal to unity and 
remains unity for greater values of ¢. Consequently, 
if the integrals are to be taken right across the bound- 
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ON POHLHAT 


and fh will be different according as A is 


To illustrate the evaluation of k, if A < 1, then A < 


b, 


l. 


ind so we must integrate with y running from 0 to 


Hence, 


. *K *] 
fy dn / fg dn + | fg dn = 
J! 0 JK 


<_ 


me 
/ fg dy + | / dy 
JO oF K 


since g = lfor¢> 1 i.e., for n> A. 
On the other hand, if A > 1, then A > 6, and ¢ must 


run from 0 to | so that 7 goes from Oto A; hence, 


"A re “A 
h = / gdn / fg dn + / fg dn = 
( 0 1 


. . 


7 
| |g dy 
Pe A I 


since f = Ofor n > 1. 


The other integrals are evaluated in the same way. 
We obtain, for A < 1 (A > Lis found not to be ap- 


plicable here), 


g¢ = 


( = 


1/1,6S0, 6 — 13,360 


ISK /35, f= —@ = 1/38 


We note also that 


(0) = —'/,, 2° (0) = 1! 


REDUCTION OF THE EQUATIONS 


Eq. (7) simplifies by writing / = \* r', and we ob- 


tain 


JO 


O'S; 


Fic 


] 














02 OF O4 O5 06 O7 O08 O9 10 
. 
Boundary-layer thickness 6; = (Q/ Ro? sin a)'/*5. Present 
paper Taylor - - - -- - 


SEN’S METHOD 180) 


i 18(35 — 72K)/d5r = (336 — 691.2K)/r 
= (a — BK)/r (say) (9 


Eq. (S) may be written 


., —10,080 — A(4I + 3/’r) 
Kk” = (10 
120r/B 


where 


A = K*(9K? — 35K + 42) | 
B = K2(12K? — 35K + 28)} 


after some more fairly heavy algebra. 

We integrate these equations with the initial condi 
tion / = Owhenr = 1. 

Eq. (10) has a singular point at the start, since the 
denominator on the right-hand side vanishes. Since 
physically we expect A’ to be finite, we must have 


10,080 + A(4/ + 3I’r) = O 
Hence, at r = | using Eqs. (9) and (11), we obtain 
K(72K — 35) (9K? — 35K + 42) — 350 = 0 


for the initial value of A. 

This gives A = 0.87933, which is the only root lying 
between Oand |. 

It is also necessary to find the initial value of A’. 
To do this we use Eq. (10), whose right-hand side has 
the limiting form 0/0; so, to find the limit, we differ 
entiate numerator and denominator with respect to 7 
and then let 7 tend to unity. 

We obtain 


K'l14 1 dA 4 ISA | A 
. 0B dK ' 10B(72K —35)} —-35B 
which gives A’ = —0.06824 at the start. 


THE SOLUTION OF THE EQUATIONS 


K changes slowly, and so we use a step-by-step 
method for Eq. (10), assuming A is constant through- 
out each step. In Eq. (9), because of the large multi- 
plier of A, this is not sufficiently accurate, but we can 
improve the accuracy by assuming throughout the 
step 


K — Ko = Ko'(r — fr) 


where the suffix 0 refers to values at the beginning of 
the step. 
Integration then gives 


I, — Io = (BKo — a — BKo'r) (log ro — log r;) + 
BKo'(ro — nr) 


where a and 8 have the values given in Eq. (9) and 
where the suffix | refers to values at the end of the step. 
Improvement at every other step was made by the 
method of Duncan.'® The steps were made at inter- 
vals of 0.1 for r from r = 1 tor = 0.4, and thence by 
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steps of 0.05 down to r = 0.2. It appears probable 


that three-figure accuracy is obtained. 


The results are shown in Table 2 and Fig. 1, in which 
Taylor’s results are shown for comparison. It will be 
seen that the thickness of the boundary layer de- 
creases much more quickly than that of Taylor as we 
approach the apex. In Table 2 and Fig. 1, we have 
written 


6; = (Q/R,? sin a)'6 


We show also values of x, the angle the surface fric- 
tion makes with the generators of the cone. 

It is of interest to note that, in the solution given 
here, 


quite closely. 


DISCUSSION 


The main purpose of this paper was to show that it 
may be incorrect to take equal 6’s and to give an illus- 
tration of a problem that has now been solved both 
ways. 

It is not, of course, certain which of these results 
more closely approximates the exact solution. Pohl- 
hausen’s method in the form applied here has been 
shown to give accurate results in other cases of non- 
retarded flow; therefore it seems likely that the result 
given here is more accurate. The algebraic and arith- 
metic labor is greatly increased however, and some dif- 
ferent approximate method must be found if it is to 
have any wide application. A more serious consider- 
ation, however, is that so far no solution of any three- 
dimensional boundary-layer problem has yet been given 
which has not some form of symmetry somewhere. 
All solutions known, exact or approximate, have been 
such that a coordinate system could be found in which 
all derivatives with respect to one of the coordinates 
vanish. One exception is that of Carrier,’ but that 
paper is believed not to give the solution of the problem 
it sets out to solve. 


The problem of edge effects, of which Carrier’s prob- 
lem is a type of illustration, also urgently needs con- 
sideration, and so far only one attempt at this has been 
made (by Howarth'*) which gives the results of great 
qualitative interest for the uniform flow part on in- 
finite quarter plate. 


TABLE 2 


, 20 0.95 0.9 0.8 0.7 0.6 0.5 
K 0.879 0.883 0.886 0.893 0.902 0.912 on 
I 0 14.44 28.9 62.7 100 144 199 
5 O 1.85 2.09 2.25 2.22 2.08 1.8 
x 90° 60°8 . 51°6 10°5 34°91 =28°7 24° 
y 0.4 0.35 0.3 0.25 0.2 

K 0.940 0.949 0.962 0.975 0.994 

I 267 309 359 $2] 197 

6 1.62 1.47 1.31 1.138 0.94 

x 213 19°8 18°2 16°7 15°32 
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RIEF REPORTS of investigations in the aeronautical sciences and discussions of papers pub- 


lished in the JOURNAL are presented in this special department 
as quickly as possible after receipt of the material. 
sponsible for the opinions expressed by the correspondents 


words in length 


The Stability of the Laminar Boundary Layer 
in a Compressible Fluid for the Case of 
Three-Dimensional Disturbances* 


D. W. Dunn and C. C. Lin 
Massachusetts Institute of Technology, Cambridge, Mass. 


March 21, 1952 


N THE STABILITY THEORY for an incompressible fluid, Squire 
| has shown that the stability problem for three-dimensional 
disturbances is equivalent to that for two-dimensional disturb 
ances. If a similar attempt is made for the boundary layer in a 
compressible fluid by writing down the linearized equations for 
three-dimensional disturbances with the type of approximation 
described in reference 2, it is found that the final equations will 
transform to corresponding two-dimensional ones except for three 
terms in the energy equation.t| However, a more careful ex 
amination of the orders of magnitude of the disturbances them- 
selves shows that many of the terms retained in the previous 
method, including the three troublesome ones, are just as small 
is those discarded, at least for subsonic free-stream velocities 
and possibly for ordinary supersonic velocities. The further 
simplification of the equations obtained by dropping these terms 
permits the establishing of a generalization of Squire’s theorem 
rhe present account is a preliminary report of this result 

The final, approximate, dimensionless, and linearized equa 
tions for the amplitudes of small disturbances are as follows 


ap |i(% of + we (ia/yM?)r + (mi/R)f" (1 
a’*p|i(w o¢g| = (1 7M? T (2 
. ( Je ) ” » 
ap |i(w c)hh| = (i8/yM*)r + (m/R)h (3) 
aly c)r + p'ag pl la + 78h + ag’) = 0 (4) 
vp |i c)0 T’¢]| = (¥ L)oT(iaf + iBh + ag’ 
(¥ Ro 0" (5) 
r/p = (r/p) + (0/T) (6) 


The notation is the same as in reference 2 except for the quan 
tities h and 8 corresponding to the three-dimensional character 
of the disturbances. The appropriate solutions of the exponen 


tial type are now illustrated, for example, by 
u;’ = h(x2) exp [i(ax, + Bx acl 
for the third component of the velocity disturbance. 


These equations are correct except for terms of order (aR 


at least for 17 <1. At higher Mach Numbers, some of the re 


*One of u D. W. D.) would like to thank the Research Council of 
Ontario for the financial assistance of a scholarship during the period of this 
research 

t This difficulty has also been noticed by others—e.g., Profs. L. Lees and 
M. Lesser 


Publication will be completed 
The Editorial Committee does not hold itself re 
Contributions should not exceed 800 


jected terms may become important. However, tentative in 
vestigations suggest that Eqs. (1) to (6) may still be valid up to 
V7 = 2, say. As may be expected, it is easily verified that the 
approximate solutions of reference 2 are solutions of the two 
dimensional versions of the present equations to the same degree 
of approximation. A restriction on the validity of these solutions 
has been pointed out in reference 3. 

Clearly, Eqs. (1) to (6) transform to corresponding two- 


dimensional ones by the transformations 


af = af + Bh am = am 
ag = a¢ C = < | 

ant = ar a =@vWy a® + B? (7) 
a? = a M = (a/a)M \ 

ad = ad R = (a/a)R 


This is the formal generalization of Squire’s theorem to the com 
pressgble Case However, since the distribution functions of the 
mean quantities depe:d on the Mach Numbers (and possibly also 
on the Reynolds Number), the physical interpretation is not 
simple 

The case of high Mach Numbers is more difficult, mainly be- 
cause of the rapid variation of some of the mean quantities with 


M. This is shown in the case of 7, for example, by the following 


result for a boundary layer with insulated wall and Prandtl 


Number of unity: 


1)/2|AM*(1 u? (8 


Tes 
Evidently, the correct order of magnitude will be given by this 
for any boundary-layer flow with Prandtl Number and wall tem 
perature not too different from their values in this special situa 
tion. It is clear from Eq. (8) that, for large enough Mach Num 
ber, 7 can become extremely large and change by a whole ordet 
of magnitude across the boundary layer. This makes it difficult 
to estimate the orders of magnitude of both the mean quantities 
and the disturbances for large JJ. Present indications are that, 
for higher and higher Mach Numbers, more and more terms pre 
viously neglected become important and that, finally, the equa 
tions become so complicated that the transformation to the two 
dimensional case becomes impossible. It should be noted that 
for such high Mach Numbers even the boundary-layer theory it 
self may require modifications 
Further investigations are being carried out to determine the 
appropriate equations and the upper limit of validity of the 


above transformations in the supersonic range 
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Some Considerations of the Stability 
of Laminar Parallel Flows 


Martin Lessen 
The Pennsylvania State College, State College, Pa. 
March 5, 1952 


’ | ‘HE TWO-DIMENSIONAL DISTURBANCE EQUATIONS for a plane 


quasi-parallel laminar flow are: 


Continuity: 
ia(U — c)p + pliau + v’) + p'v = 0 
Force: 
plia(U — chu + U'v) = —iap + 
rs ” 4 ° ta , , 
Sin = a8 + — 1 + plu + la 
3 3 
: : ie. : la , 
alia(U — cv] = —p + al -0" —av+—4 | 
3 3 
4 2 
a’ y’ f iaw ) 
3 3 
Energy: 


peylia(U — c)T + T'v] = R(T" a?T) + K’'T’ 


Pliau + v') + QU "(wu + sai 
where 

U(y) + u(y)eta*— © = velocity in x direction 

v(y)etat® —) = velocity in y direction (1 to 4 
direction ) 

pt ply’ —? = pressure 

T(y) + T(y)ete*— = temperature 

p(y) + ply \e tals — cl) = density 

By) + w(yete*— = viscosity (absolute) 

K(y) + K(y)e’* ct) = thermal conductivity 

a = disturbance wave number 

r = disturbance velocity of propaga 
tion 

l = time variable 

C = specific heat at constant volume 


(for a gas) 
Consider first the case of an incompressible fluid. For this 
case, the energy equation need not be considered. The equations 


become 


Continuity: 


iau +v' = 0 
Force: 
pliaelU — chu + U'v] = —-iap + pu” — a®u) + p’(u’ + av) 
plia(U — cv] = pb’ + wv” a’*y) a’ (Zien) 


For the problem of incompressible flow with variable viscosity, 
if p is eliminated between the force equations and wu eliminated 
by continuity, the resulting equation in dimensionless form is, 


(l -chv — av) U"vy = (7 aR) {av - 2a*v T a’v) + 


The boundary conditions for the foregoing equation are that 
v and v’ vanish at the boundaries. 

Consider now the case of laminar ‘‘mixing’’ of dissimilar fluids 
in which the properties of each fluid remain constant 

The density and viscosity are step functions in y, the dis 
continuities occurring at the interface between the two fluids 


The equations therefore become 


Continuity: 


tau + 1 = (0) 


Force: 
plia( clu + U'y| = lap + plu” an 
plia( U CN = p’ + pv” av) 


Reducing the above to a single equation in v as before, th 


resulting dimensionless equation is 


" I\ 
U"ys = (1/aR)r(i 2 


oat + a 


where the characteristic density and viscosity are taken as 


and j respectively, and v is [(%/m)(j1/p)]. The boundary 


conditions on v are v = v’ = O at the boundaries 
It is necessary also to prescribe four conditions at the interfac, 


between the fluids in order to integrate the disturbance throug} 


it. These conditions are 


where the subscripts pertain to the contacting fluids. Conditio 
c) is tantamount to 


f(y’ + fav,) = polite’ + fave) 
or 
14" = (1) / fei . an {1 ( (41 / 2) 


Condition (d) can be reduced from the first force equation 


pi lial l cu, + Unt] ilu,” au,| = 
polia(l C)ts. + Uje've jo |u a 
or 
M f # p2 P 
= i a ') a(l ( 
ii jy 7D 
la 
(Un'a — 
jie 


where l; is the steady velocity at the interface 

In considering a compressible flow, it will be assumed that th 
inviscid disturbance motion is isochoric. This is tantamount to 
assuming that an element of the fluid does not change volume 
or that the pressure perturbation is of higher order than tho 


velocity perturbation. Continuity therefore becomes 


lau + 7 = UV 


la( pu) + (pi pv = U0 


If pressure is eliminated between the force equations 


continuity introduced, the resulting disturbance equation 1s 


l C)\ Cpt (pt ’ a*( pi (pl \’2 = 


The boundary conditions are again that v and v’ vanish at th 
boundaries 
a result of the Lees-Lin work—-namely, cy = UL’ where (pl 


If the assumption of an isochoric disturbance is made, Squire s 


s , +} 
theorem can be extended to the compressible case directly with 


out considering the energy equation 
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It should be noted that agreement is obtained with 
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Taking Account of the Compressibility of the 
Material in the Plastic Buckling of Plates 


pD Pp. Bijlaard 
Professor of Structural Engineering, Cornell University, Ithaca 
January 15, 1952 


. A READER'S FoRUM ITEM in the JOURNAL OF THE AERONAU 
rICAL SCIENCES of November, 1951, Stowell and Pride give 


formulas for plastic buckling of plates which are said to take ac- 
count of the compressibility of the material and compare the re- 
sults with those from Stowell’s earlier computations.' How 

ever, as was explained in reference 2, Stowell’s previous formulas 
follow directly from the present writer’s earlier theory,** from 
which they were derived via Ilyushin by means of the simplifying 
issumption that Poisson’s ratio is equal to 0.5—that is, the com 
pressibility of the material was neglected at that time. Hence 
Stowell and Pride should now arrive again at the present writer’s 
original formulas. It follows directly from their Fig. 1 that 
they do not, since for both cases they find a correction with 
respect to Stowell’s earlier results which is opposite to that follow 

ing from the writer’s theory. 

This leaves the question of which formulas are right, those of 
Stowell and Pride or those by the writer. The former do not 
give any derivations, while the writer derived his formulas in 
several different independent ways.**° It appears that Stowell 
and Pride’s formulas may not be correct in view of tie following 
irguments 

Consider first the buckling of an infinitely long compressed 
flange, simply supported at one of the unloaded sides The 


stress-strain relation 
[21 + uw (1 


where 
= (1/2 .1(1/2 AE E! 9 
M ( « r {4 “ L User Li 


is used by Stowell and Pride, was derived in reference 6. Here, 
vis Poisson's ratio for the elastic deformations. With incipient 
buckling, only infinitely small shear stresses are superimposed on 
the initial purely compressive stresses. According to the theory of 
plastic deformation, as used in all pertinent derivations, the re 
sulting shear strains are the same as if the compressive stress and 
the infinitely small shear stresses would increase in the same ratio 
is that of their final values. The infinitely small shear stresses 
do not increase the equivalent stress, and, hence, the secant modu 
lus is the same as that just before buckling, as denoted by Ege, 

Hence, the infinitely small buckling shear stresses are, from Eq 


1), 
dr = (dy)Esec./[2.1 + w)| = EF dy (3 
instead of dr = Gdy in the elastic domain. In the elastic domain, 


the buckling stress is 


where the plate flexural rigidity 

N = Ef?/(12(1 ve 
Hence, in the plastic range the buckling stress is 
Crp = (EF/G)o-ry, = } rk Ege, {[12(1 + w) (1 v) | f(t b)? (4 


This is (1 a)/(1 v) times the expression given by Stowell 


ind Pride. The reduction coefficient for plasticity is 


7 = EF/G = {(1 + »)/(1 4 pb) ( Esec./E) (a) 


This is in accordance with the present writer's earlier results." 
Since w is larger than vy, 7 is smaller than Eg,, 
compressive curve for this case falls below the stress-strain dia 


gram, instead of above it, as in the note by Stowell and Pride 


E, so that the plate 
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For the second case of a compressed plate, simply supported 
at the unloaded edges, the simplest way to see which formulas 
are right is to assume a stress-strain diagram as in Fig. 1 and to 
consider buckling at point A, where the buckling stress o; (Fig 
2) is equal to the yield stress ays, Eyer. = E and the tangent modu- 


lus EF, = 0. 


8 


With the relations *: 
doy = E(A der + B dey) | 
do, = E(B de, + D de,) > (6 
drry = EF dyz \ 


between excess stresses and excess strains during buckling, 
the differential equation for this case is*** 
EI [Ad‘w/dx* + 2(B + 2F) d4w/dx2dy? + D d4w/dy*] 4 

tocrp O?w/Ox? = 0 (7 


With 
w = wy cos (prx/a) cos (gry/b) (8 


and after choosing p and g such as to make o¢rp minimum, it fol- 


lows 


2r°*EI/b*t) (VAD + B + 2F) (9) 


Ccrp = 
According to the deformation theory (Hencky—body), the ex- 
cess stresses and excess strains are simply the difference between 
the total stresses and strains after incipient buckling and those 
just before buckling, all calculated as if stresses and strains in 
creased proportional to their respective final values.* Denoting 
the secant modulus for the plastic deformations alone as Ep and 
recalling that for the plastic deformations Poisson’s ratio is 
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0.5, the plastic strain in the }-direction is 
éyp = (y/Ep) — loz/(2Ep)] (10) 


while the plastic strain in Y-direction is 


€rp = (or/Ep) — [oy/(2Ep)] = [(2er ay) /(2a, 


or) le€yp (11) 


On the other hand, the plastic strain is the difference between 


total and elastic strain, so that 


€yp = €y éye = ey — [(o,/E) (vo;/E)| (12 


8 


Hence, from Eqs. (11) and (12), 


€r = €re + Exp = (o7r/E) (vo,/E) + 
((2oz — oy)/(2oy — or)) le (o,/E) + (voz/E)| (18 


Assume that the stresses and strains in Eq. (13) are those after 
incipient buckling. The stresses and strains just before buckling 
are 

Sr = ys; ty, = 0 

éx = oys/E; éy = —voys/E (14) 
so that the stresses and strains after incipient buckling can be 


written as 


Or = dys + doz, oy, = doy 


ér = (Gy E) + dey, € = (voy, k t de, (15) 
Since according to Fig. 1 the stresses have to remain at the yield 
stress, they have to satisfy the condition for constant distor 


tional energy 


o;" TrOy + Gy? = Cys’ (16) 
or from differentiating 
doy, = |(20, Oy)/\ Or 2ey )|do,z 
from which, from Eq. (14), 
doy, = 2do, (17) 


Insertion of Eqs. (15) and (17) in Eq. (13) yields® 


doy = |E/(5 tv)] (der + 2Qde,) | 


doy = |E/(5 — 4v)| (2der + 4de,)f sa 
so that, from comparison with Eqs. (6), 

A2z1i1/(5 ly), B=w2A, Dz 4A (19) 
Since for the present case Exec. = E and pw = », from Eq. (3) 
F = 1/(2(1 + »v)!] (20 

Insertion of Eqs. (19) and (20) in Eq. (9) gives 
(Gerp)B = 118 (1 + v) (5 ty) |} (7°21 /b*t (21) 
The same result follows from the general equations, Eqs. (22) and 


(23) of reference 8 or Eqs. (21) of reference 2, for 8 = 0, e = Vor 
E,; = Eandtan ¢ = Vor & = 0. On the other hand, according 
to the formula given by Stowell and Pride, for the present case, 


with Esec. = E, wp = vand E; = 0 
(ccrp)s-p = [2/(1 v)| (7? EI /b*t) (22) 
With » = 0.5, Stowell’s result in reference 1 follows from Eq. 
(21) as 
(ocrp)s = 42° EI /b*t (23) 


In the elastic domain 


Ocre = [4/(1 v?)| (rw? EI /b?t) (24) 


By comparing Eq. (21) with Eq. (24), the reduction coeflicient for 
plasticity is, with »y = 0.3, 

np = 0.83 (25) 
while from comparison of Eqs. (22) and (24), according to Stowell 
and Pride, 


ns—-p = 0.65 (26) 


rICAL SCIENCES Toes, FOSS 


On the other hand, Stowell’s earlier results follow from eo, 


parison of Eqs. (23) and (24), by insertion of vy = 0 5. Vielding 
5 


ns = V.i0 = 


Hence Stowell’s value ns has to be corrected upward to ().88 nA 
not downward to 0.66, as would follow from the formulg “ 
Stowell and Pride : 

To settle the above question it would be appreciated if Stoy, 


and Pride would show how they derived their formulas 
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_— AUTHORS AGREE that Curve 2 should be under Curve 
as pointed out by Professor Bijlaard The & in the flange 


buckling formula is a function of Poisson’s ratio and so is itself 


slowly varying function of stress. When Curve 2 is corrected 


for this effect, it is displaced downward a minute amount so that 


t 


it lies slightly below Curve 1 rhe corresponding & in the plat 
buckling formula is independent of Poisson's ratio, being equa 
to 4, and so curves 3 and 4 are unchanged. In the example cho 
sen, our calculations show oys to be the buckling stress in agret 
ment with Professor Bijlaard’s valuc 

rhe object of the note was to give aircraft designers a numeri 
whet 


rh 


estimate of the errors due to neglect of compressibility 
computing buckling stresses by the method of reference | 
obvious conclusion is that such errors are quite negligible com 
pared with the usual uncertainties due to variation in mater! 
properties 


1 


Derivations were omitted because they parallel the origi 


minut 


derivations of reference 1 and because the resulting 


differences did not justify a lengthy paper 


A Generalized Ground Run Chart 


William F. Savage 

Assistant Professor of Aeronautical Engineering, Univer 
Kentucky, Lexington, Ky 

April 3, 1952 


SUMMARY 


\ chart is presented which can be used to determine the ground rut 


tance in take-off for any aircraft having a power plant that develops co! 
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he chart is 


nt thrust throughout the range of ground run velocities 
stan 
plicable for use in preliminary performance estimates on turbo 


especially ap 
t and rocket propelied aircraft 

te 

The equation used in the development of the chart is based on the simple 


acceleration determination of 


equations normally used in the 


force and 
introduced to facilitate 


Simplifying parameters are 


round run distance 


jing of the chart 


ead 
SYMBOLS 
iircraft acceleration, ft. per sec. per sec 
iireraft drag coefficient with ground effect (all wheels on ground 
aircraft lift coefficient with ground effect (all wheels on ground 
vl maximum lift coefficient of aircraft with ground effect 
D drag of aircraft, lbs 
acceleration of gravity, 52.2 ft. per sec. per sec 
K ul Cb 
lift of aircraft, lbs 
take-off wing loading, lbs. per sq.ft 
S wing area, sq.ft 
§ ground run distance, ft 
take-off thrust of aircraft power plant or plants, Ibs 
ratio of aircraft thrust to aircraft weight 
velocity of aircraft, ft. per sec 
| take-off velocity of aircraft, ft. per sec 
WW take-off gross weight of aircraft, Ibs 


coefficient of rolling friction 
mass density of air, slugs per cu.tt 


ratio of mass density of air at altitude to that at standard sea 


level conditions 


INTRODUCTION 


I’ PRELIMINARY PERFORMANCE ANALYSIS of propeller-driven 
ground 


This method is usually found to be less time con 


iircraft, run distance is determined by graphical 


integration 
suming than direct integration, since the thrust developed in the 
ground run 1s When that 


develop constant thrust during the ground run are 


function of the velocity aircraft 
considered 
rocket-propelled aircraft), direct integration is 


turbojet and 


feasible, and a chart can be constructed which is applicable to 
| such aircraft having constant-thrust power plants 

As in the case of modern aircraft having tricycle landing gear, 
run is considered with all three wheels on the ground 
At that 
point, a nose-up attitude of the aircraft is attained by the pilot 
The take-off speed is based 


ittaining a lift coefficient equal to 90 per cent of 


ground 


until the ‘‘unsticking’’ or take-off speed V,, is reached 


ind the aircraft leaves the ground 
on the aircraft 
the maximum value. 


THE Basic GROUND RUN EQUATION 
rhe basic integral for the determination of ground run dis 
tance can be written in the following mannet 
el 
bY, { V dV/a (1 
e t) 


Considering the sum of the horizontal forces acting on the ait 


ft and equating this sum to the product of the mass of the air 


raft and its acceleration, we will obtain 
/ u( / D (Hu 2 
ol 
( P/u im ul /V (D/i 3 
Replacing the acceleration term in Eq. (1) by its equivalent 
from Eq. (3) and substituting the equivalent equations for lift 


nd dr ig, 


. 


‘ V dV 
; ; j 
/ iz. 7 (ee moe) 
2W 


iH 


Upon integration of Eq. (4), an equation will be obtained in 


terms of V,, but V,, may be eliminated by the substitution 


Vu! 2NV/0.9CiwpS 


= ; . on 
«he equation resulting from integration and simplification will be 


RS’ 
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GROUND RUN- 1000 FT. 


Fic. 1 
TABLE | 
Runway Surface m 
Concrete 0.02 
Macadam 0.08 
Short Grass 0.05 
Long Grass 0.10 
Sand 0.30 
W uCy Cp . 
Se = > . ance s + | (5 
Sgp( uC Cp) O0.9Cry (7 /U pw! 
Letting A = uC Cp, m = W/S, and t = T/W and sub 


stituting the values for gravitational acceleration and the altitude 
density ratio, the final equation will be obtained, 
K 


F Mm 
Sg = 13.04 - In . + | (6 
ak O.9C, y(t mM 


It will be noted from Eq. (6) that a short ground run distance 


can be low wing loading, high maximum lift co 
efficients and high ratios of thrust to airplane gross weight rhe 


at high thrust 


attained by a 


effect of A on the ground run distance is slight 
weight ratios and at high lift coefficients. Consequently, an ait 
plane capable of developing high lift coefficients (low take-off 
speeds 9 high gross 
at take-off, is essentially independent of the lift and drag 


and with a thrust compared to airplanc 


weight 


coefficients obtained while the aircraft is in the three-point atti 
tuck When A 0, Eq. (6) will reduce to 
m l ‘ 
Sq = 13.04 . (7 
¢ L0.9 Cim(t Me 
Eq. (7) will give minimum ground run distance (considering all 


other variables remain fixed ) 


THE GROUND RUN CHART 


Fig. 1 is the generalized ground run chart that is based on Eq 


6 Table 1 lists the values of the coefficient of rolling friction 





used in the construction of the chart. It will be noted that the 
chart is based on the no-headwind condition in take-off. Head 
wind has not been considered, since preliminary performance 
estimates are usually based on still-air conditions. 

As an example of the use of the chart, consider a turbojet 
propelled airplane having the following characteristics: take- 
off thrust, 5,000 Ibs.; gross weight of airplane, 20,000 Ibs.; wing 
area, 400 sq.ft.; Cray, 2.0; Cr, 1.0; Cp, 0.070; pu, 0.02 (concrete); 
and density altitude, sea level. 

Thus, A = —0.05, ¢ = 0.25, and m = 50. 

Entering the chart with the specified value of ¢ and proceeding 
as shown by the dotted !ines, a ground run distance of 1,680 ft 
is obtained. 

It is believed that the use of the generalized ground chart will 
greatly reduce the time required in preliminary performance anal- 
ysis. It must be emphasized that, of necessity, the chart shown 
herein is greatly reduced in size, but great accuracy can be ob 
tained with a large chart 


Note on the Drag of a Finite Wedge at Mach 
Number 1 







Isao Imai 
Department of Physics, Faculty of Science, University of Tokyo, 
Tokyo, Japan 


March 20, 1952 


iat or bore & 


prea GUDERLEY AND YOSHIHARA! made an elaborate 
study on the flow over a wedge profile at Mach Number 1, 
using the hodograph method as simplified by the transonic law of 
similarity. With the same approximation, Cole? investigated 
the velocity distribution on the front part of the wedge at high 


Pipa tee 


subsonic speeds and obtained analytic expressions for the nose 
drag of the wedge. 


sir 
girr ei 


In this note it will be shown that a simple consideration using 
the idea of Kirchhoff’s dead-air theory gives a satisfactory result 


t= 


for the nose drag of a finite wedge. 

The wedge is assumed to be symmetrical about the undisturbed 
uniform flow and to have the semiangle 6 (Fig. 1). The stream 
lines flowing along the wedge are assumed to separate at the cor- 


‘fe aie 


ners A and A’, making free boundaries AB and A’B’. 
The exact hodograph equation can be written in the form® 


‘ 


See ia 


(O2y/dr2) + (02y/062) = R(r)y¥ (1) 
where y = A'/‘W, K = (upo/p)?, »w = (1 — q? cyt 7 = 
q a 9 , ° 

-f, ug! dq, and k = K~'/‘d?K'/‘/dr?. Here, q and @ are 

the magnitude and direction of velocity vector, p is the density, 

po is the stagnation density, ¢ is the local sound velocity, and W is 

the stream function, the critical velocity being taken as the unit 

of velocity (qx = c* = 1). For adiabatic flow, if the velocity is 
close to that of sound, we have, approximately, 


(4 











RB’ 


Fic. 1. 
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k(r) = —(5/36)r 
where y is the ratio of the specific heats. With this “transg 


approximation,” it can readily be seen that, if Y = yo(r, 8) is 


solution of Eq. (1), then y = yo(ar, a0), where a is an arbitrary 


constant, is also a solution of Eq. (1). (This forms essential 
the basis of von Karman’s transonic similarity law 
It can be shown without difficulty that the discontinuous fio 


past a flat plate accompanied with a dead-air region is, 


in gen 
eral, represented by 
: QV» (q) : 
Get F «a sin 2n@ 
—os 
l Qeilg 
where C is a certain constant and Q,(q¢) = K W (7), wher 
W,(7) is the solution of 
(d?W’,, /dr*) n*W, = RW 
such that W,(7) ~e>"7 as r= 
With the transonic approximation k(7) 5/36 )r h 
relevant solution is 
Q0.(q) = K-"'/*7'/2H1,,(nir 
Hi, being Hankel’s function of order ! For simplicity, 


only the case of Mach Number 1 will be considered her Thus 


To — 0, so that 
Hi, (nitro) ~ i wr 41/3) (nr2/2 
and Eq. (8) becomes 
y= K’/W=(C, > (in)'/87'*Hiy,)(2nir) sin 2n0 { 


C, being another constant related to C. From this expression, 
the stream function for the flow past a wedge of semiangle 4 
can be immediately obtained by the transformation 


er 


r—>7rr/26, 0 10/256 
Thus, 
y= C. > (in) “7 ?H “\(inar/6) sin (nx0/6 5 
C, being another constant. 
The line element of the stream line ¥ = 0 is, in general, giver 
by’ 
ds = g°'K'/*|(Oy/00) dr (Oy /Or) dé 6 
which becomes, for the wedge profile 6 = 4, 


ds = C y —1)0n' 97? /87/9TT), Ciner/b) dr 


For vanishingly small semiangle 6 —~ 0, Eq. (7) reduces, because ol 
1 


the asymptotic expansion for the Hankel function // , to 
ds =€,>,( 1)" ‘% se as dr ) 


C,; being still another constant 


Now the pressure coefficient is defined as 


: p p 2 p 
(1/2 Pod« Y p 


with 
p/po = [1 + (4 1) (1 q*)/2 
Hence, with transonic approximation g = 1, 
Cp (1 q*) (y + 1)/2 (37 by 


The drag coefficient of the wedge is obviously given by 
Cp = sf Cy ds J ds 10 


Therefore, remembering that the nose and corner of the wedge 


correspond, respectively, to 7 = © and rt = 0, we have 
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F(11/6, 2/3) 
where 
° 
Fla, b) = { p (—1)" n2 e—977/8,0—1! dy 
J 0 1 
= (6/r)’T(b) ¥ ( 1)"n2—° 
= (6/r)’R(b) (22 o+t 1) ¢(b a) 


where ¢ is Riemann’s Zeta function. The successive transforma 


> ( and b > a, but the final result is 


tions are justified only if b 
valid without the restriction b > a because of analytic continu 


ition. Finally, after some simple transformation, the reduced 


drag coeflicient of the wedge is obtained in the form 


Co = (7 7 4) "8 Cp 
2°93 2 1 ¢(3/2) 
= = = ———— oor = 1.488 
7 2°/7* — 1 $(13/6) 12 


Somewhat similar calculation based on the simple approxima 
tion k(r) = 0, which has been adopted in the author’s previous 
paper,’ gives 

2/ - | (7 6) T 


Co = 2 X 3° — tan = 
i 1 ¢(11/6) 12 


1.944 


These two values are to be compared with Guderley-Yoshihara’s 
numerical result 


Cogy = 1.75 


1 — 2/3 2 
(r(2/3)}2° (7) - 


This is probably due to 


and €ole’s analytical result 


toe = 4(5) 
Cpe = 
. 3 


The agreement is surprisingly good. 
the fact that the flow condition downstream of the corners of a 
finite wedge has only a little influence on the pressure over its 
“‘dead-air”’ theory 


1.67 


front part. This seems to indicate that the 
might be useful for the consideration of nose drag of bodies at 


high subsonic speeds. 
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Note on the Conditions at a Sharp Leading 
Edge in a Supersonic Stream* 


Russell E. Duff ¢ 
Department of Physics, University of Michigan, Ann Arbor, Mich. 
March 20, 1959 


I A RECENT ARTICLE, Bardsley! has described the flow con- 

figuration in the vicinity of the sharp leading edge of an 
airfoil at an angle of attack of 15° at a Mach Number of 1.965 
The configuration observed is shown in Fig. 1. In the region 
* The author wishes to express his gratitude to Robert N. Hollyer, Jr., 
lor his active assistance in this research program and to the Office of Naval 
Research for financial support furnished through contract N6-ONR-232 TO 
IV 

t Now at Los Alamos Scientific Laboratory, Los Alamos, N.M 
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above the wedge where inviscid theory predicts only a Prandtl 
Meyer expansion wave, a weak shock was found preceding this 
wave. Bardsley explained the phenomenon by considering 
only the bluntness of the wedge 


ration and boundary-layer growth did not contribute to the for- 


It was shown that flow sepa- 


mation of the leading shock under these conditions 


However, recent shock-tube experiments suggest that this 
configuration is by no means unique for a sharp airfoil at an 
angle of attack in a supersonic stream; rather it seems to be the 
At high subsonic Mach 
of attack 


that is, 


limiting case for high Mach Numbers 
Numbers, the upper surface of an airfoil at an angle 
is covered by a detached, turbulent boundary layer 
the airfoil is in a state of stall. At higher Mach Numbers this 
boundary layer becomes thinner, and at some transonic value 
the flow reattaches to the rear portion of the airfoil. The region 
of separated flow shrinks rapidly as the Mach Number is raised 
further and presumably ceases to exist at some value below 
1.965. 

Fig. 2 is a shock-tube shadowgraph of the transonic flow in 
nitrogen about a double-wedge airfoil at a Mach Number of 
0.85. In this case the length of the separated boundary layer, 
or the compressibility bubble, is approximately '), in 

A photograph showing the flow configuration at a 
Mach Number is presented in Fig. 3. This picture is a shadow- 


higher 


graph of the flow in carbon-tetrachloride vapor at a Mach 
Number of 1.65 about a sharp 12° wedge at an angle of attack 


Carbon-tetrachloride vapor was used in order to obtain 
Mach 


of 15°. 


relatively high Number without sacrificing resolution. 










WEAK SHOCK WAVE —~, 
PRANDTL- MEYER 
WAVE 
SHOCK WAVE 
Fic. 1. Supersonic flow at the tip of a wedge, as observed by 


Bardsley. 





Transonic flow 


Fic. 2 showing a large compressibility bubble. 
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Supersonic flow observed in a shock tube showing weak 
shecks at both sides of a Prandtl-Meyer wave. 


Fre. 3. 


The effect of temperature on the specific heat of the vapor is not 
important in this case, and elementary shock theory is adequate 
for describing the flow in the shock tube. 

Fig. 3 shows a shock wave at both the front and rear of the 
Prandtl-Meyer wave. The first shock is caused by the bluntness 
of the wedge as described by Bardsley. 
by the compressibility bubble as shown in Fig. 2. However, 
the bubble itself is too small to be seen distinctly in this picture 


The second is caused 


It is interesting to notice that particle paths can be seen below 
the wedge in Fig. 3. Apparently, the disturbances caused by 
small irregularities on the windows of the tube are made visible 
by the very high index of refraction of carbon-tetrachloride 
Similar paths are not found when other gases are used 


vapor 
in the shock tube 
Clearly, the wake in Fig. 3 is still in a transient state. The 


role of several of the observed shock waves in the establishment 
of a stationary flow from the unusual initial conditions en 
countered in a shock tube will be discussed in a forthcoming 


paper 
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Effect of Flapping Hinge Offset on Rotor Blade 
Stall 


J. P. Chawla 
Hughes Aircraft Company, Culver City, Calif 
March 24, 195¢ 


INTRODUCTION 


M" LER HAS SHOWN in his analysis of helicopter control and 
stability in hovering flight! that stability and control may 
be improved by incorporating an offset flapping hinge in the 
rotor design. This improvement is due to the fact that the aero 
dynamic and inertia shears act at a greater arm when the flap 
ping hinge is located well out along the blade radius, thus yield 
ing improved control and stability of the machine. 

The following analysis is submitted as a contribution to the 
in particular, the effect of flapping hinge 
It is concluded from 


design of such rotors 
offset on the 
this study that, in the helicopter régime of flight and with normal 


tall at retreating blade tip. 


AERONAUTICAL 


SCIENCES Ty, LSSzZ 


values of Lock’s Inertia Number y = (pacR')/TJ,,, the stall condj 
tion at the retreating blade tip is not materially improved by 
. J 


hinge offsets up to 20 per cent of the blade radius 


THEORETICAL DEVELOPMENT 


Inertia and aerodynamic forces acting on a blade element at 
distance r from axis of rotation are shown in Fig. 1 (weight of the 
blade is neglected). The blade, acted upon by these forces, js , 
equilibrium at a flapping angle 6 = A(y). 
dition of vanishing moment about the hinge as 


We express the con 


M,+MM; =0 


or 

*| | 

| R(x e)dL = [Q?R2Bx(x — €) + R2B(x e)?] dy 
. 7 € r € 
Let 


/ 


B, cos Y — By sin ¥ Qt 


(x e)? dm 
r € 
3 
¢ = (eR? 1) | (x e) dm 
I € 


Substitution in right-hand side of Eq. (1) gives 


ll 
~ 
— | 
. 


VM; = 1,2?[Bo + £8] 


For e = ¢ = O, we get the well-known expression for inerti 


moment of a fully articulated rotor, 


MoMENT OF AERODYNAMIC FORCES 
For simplicity we assume a rotor with first harmonic flapping 


and feathering control. If we substitute 


6 = Oy 4, sin y¥ 

B = By B, cos yp Bz sin y 

dL = (1/2) paca,l rR dx j 

a = 6 (Up/U7 

Up = AQR (x 6)RB — p2QRB cosy [Ais inflow 
ratio 

Ur = QR(x w sin y¥ 


and retain only the steady state and first harmonic terms in tl 


left-hand side of Eq. (1), we get 





dL (AERODYNAMIC FORCE) 











2 
xRQdm (C.F.) 


xROBdm (NORM. COMP. 
\ OF CF) 


(x-€)RAdm (INERTIA FORCE) 
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FLAPPING HINGE 


AXIS OF ROTATION 


>be |] 
ini 


= HINGE OFFSET, FT R = BLADE RADIUS, FT. 
= FLAPPING ANGLE,RAD.; COSf=!, SINB=8 
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Forces acting on blade in equilibrium at flapping ang]! 
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O?R*[ {@(A + Bu?) + ACu + BiBeu 


Vv = (1 2 )pac 4 
ib r»C} + {2uCo — [A + (3/2)Bu?] — Bi[D 
(1/2)By?} 2ryB} sin y + } B2[D + (1/2)Byu?] 
BoCu{ cos ¥ } (5a 
where 
l l ] ; 
= € € 
A { 12 
| l a 
R= € é€ 
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é, (Sb) 
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~~ s £ 
] 2 a. l 
D= € . a 12 é! 
4 } y 2 


We define ¥ 


in Eqs. (3a 


= (pacR*/I,) and equate the coefficients of like terms 


and (5a). The following explicit expressions for the 


flapping constants 8; (7 = 0, 1,2) result: 


1 es 7" I . 
: E}14 D By?) X 
Hi+¢) al 1c? 2 
4 y* . I { 
ss 7 'E F ] 
ie = a (0 bu) (6a 
hi Al(l + ¢$) ¢ 2 


7 l j ¥" . l . 
. . - CEy | D Bu? 
% aAl4(1t+¢e 2 


yeBC ly 
f l pe? 
47(1 + ¢ | 


where 
' 7? (p l B (v l B ‘) vy7eBC 
A t 3” + Be e 
i? 2 2 1((1 + ¢) 
E =0(A + Bu?) + Cu AC 
fF 2(AoC \B)u A{/A + (3/2)By? | 
We note that substitution of e = ¢ = O in Eqs. (3a), (5a), and 
5b), before performing the above algebraic operations, yields 


the well-known expressions for the flapping constants of a fully 


uticulated rotor, referred to the rotor disc plane. They are 


y | % l 
Bo = Cl gr} 4 (Gin \) 
214 3 
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By 
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APPLICATION--EFFECT OF FLAPPING HINGE OFFSET ON BLADE 


STALI 

We apply the theory developed above to calculate the resultant 
angle of attack at the tip of the retreating blade (x = 1, ¥ 
270°) when the machine is trimmed for no longitudinal tilt of the 
rotor In other words, 6; is so chosen as to make £, = 0 
@tip, y 270° is then given by 

@)tip, ¥ 270 A +4 — [A/C M f 
where 
8/3 )ulé (3/4)A] 
A 
l (3 2)? 
lore = 0, and 
j ( Cl@(A + By?) AC \t 
a 2(A ¢ \B ( ) 
| l+¢ [D + (1/2)Byu*]  § 
wae (8 
ee ¢ C2y? 
| Bu?) 4 
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Fic. 2. Effect of flapping hinge offset on stall of the retreating 
blade tip 


Results from Eqs. (7) and (8) for a single-rotor helicopter with 


the following characteristics are plotted in Fig. 2 


Cz = (7 = W)/(rpQ?R*) = 0.0047 
A 0.0075 (f = flat plate area representing fuselage drag, 
and A = rk?) 
a = nc/rR = 0.05 (rotor solidity) 
a (dC, /da) = 6.0 (lift curve slope of rotor blad¢ 
Cp, = 0.012 (profile drag coefficient for blade 


The nondimensional parameter ¢, which is a function of the 


mass distribution along the blade, may, in general, be written as 


¢ = k[e/(1 € 


k 1.5 corresponding to a uniform mass distribution has been 
used 

For ¥ 12 (« 0), the effect of hinge offset on the coning 
ingle 8) was found to be negligible. The effect of reverse flow 


has been neglected but will be included in further studies to be 


r¢ ported later 
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Comments on ‘“‘Uniformly Loaded Semi- 
Infinite Wedge-Shaped Plates’’—Author’s 
Reply 


A. J. A. Morgan 

Aeronautical Engineering Research Consultants, Inc., Pasadena 
Calif 

April 28, 1952 


I IS DIFFICULT TO DETERMINE from Professor Reissner's re 
marks,! whether he is alluding to the simplicity with which 
the whole problem (partial differential equation plus boundary 
conditions), or the simplicity with which the ordinary differential 
equation obtained by the present author |Eq 7), reference 2}, 
can be solved. It will be shown subsequently that the latter 
certainly cannot be objected to on the grounds of ‘lack of sim 
plicity."’ The former is then reduced to the problem of deter 
mining whether the boundary conditions can be satisfied with 
greater simplicity by using Professor Reissner’s suggested solu 
tion-form (w = r‘W(@)) or by use of the solution given by Eq 


(8).2 Since, in either case, this involves direct substitutions 
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into the boundary conditions and the subsequent solution of 
four linear algebraic equations (for the general problem), it is 
clear that both methods! ? are inherently possessed with the 
same degree of complexity—the question of simplicity is then 
a matter of personal opinion and, hence, a moot point. 

The solution of Eq. (7),? can be obtained by elementary meth- 
ods, and no additional information such as“... the known fact 
that the assumption* w = r"W(@) reduces the expression V4w to 
an expression LW(@), where the operator L has constant co- 
efficients,”’ as stated by Professor Reissner, is required. An out- 
line of the method of solution for Eq. (7),? clearly illustrates this 
point and indicates that this ordinary differential equation with 
variable coefficients is inherently no more difficult to solve than 
one with constant coefficients. Let ZL denote the differential 
operator defined by Eq. (7);? then an elementary computation 
gives the following result: 


L(y") = fo(m)n"™~* + fim )n™ ~2 + fo(m)n™ (2) 
where 
fom) = mim — 1) (m — 2) (m — 3) / 
fi(m) = 2m(m — 1) (m — 3) (m — 4) > (2) 
fom) = (m — 1) (m — 2) (m — 3) (m —- 4)\ 


From Eqs. (1) and (2) it is immediately seen that » and »* must 
be solutions of the equation L(f) = 0. Furthermore, it is also 
obvious from the above that, for appropriate ratios a/b and c/d, 
the two linear combinations an! + bn? and cn? + d must also be 
solutions of L(f) = 0. The particular solution of L(f) = P/D is 
obtained by inspection. This shows that only purely elementary 
methods have been used to solve the equation L(f) = P/D. The 
simplicity of the above method is evident. 

The question as to whether or not the restriction to wedge 
angles less than 180° is necessary is still subject to detailed in 
vestigation, since the solutions obtained by different methods 
must be transformable into each other. The difficulties con 
tained in one method must then also be contained in the other. 
The use of polar coordinates alone cannot eliminate a funda 
mental difficulty. 

Professor Reissner has apparently failed to notice (by his con 
sistent use of the word ‘‘assumption”’ with respect to the solu 
tion forms w = r‘*W(@)andw = r”*4W(@)) that reference 2 serves 
as a vehicle for the introduction of a new method of approach 
to the solution of such problems. It was one of the purposes of 
reference 2 to show, by an example, that such a word can be com 
pletely omitted from presentations of this type. The required 
form of the solution, irrespective of the system of coordinates, 
can be obtained with mathematical certainty by use of group 
theoretic arguments, and known theorems* * arising therefrom, 
as indicated by Eqs. (2)-(6) of reference 2. This shows that the 
solution-forms assumed by Professor Reissner (and known to the 
author) can be classified as similarity solutions—in particular, of 
Vw = (P(x,y)/D). In group-theoretic language, such solutions 
would be classified as invariant solutions of the partial differential 
equations under consideration. This procedure* * enables one 
to consider all such problems from a unified point of view and 
obviates the necessity for knowing the required solution-form in 
advance. It has the practical advantage of enabling one to 
search for such invariant solutions in a systematic manner. This 
has also been emphasized by Birkhoff,® ® who at the time was not 
in possession of the general theorems given in references 3 and 4. 

The fact that the existence or nonexistence of such solution- 
forms can be determined (without their first being assumed) by 
the above mentioned methods* ‘' may be illustrated by use of 
reference 2. The whole argument, Eqs. (1)-(16) of reference 2, 
could have been made by omitting Eqs. (7) and (8). The sen- 
tence at the end of Eq. (16) would then read “ . . . the problem 
can then be completely solved in terms of the solution, f(7), of 
an ordinary differential equation."’ That such an ordinary dif 


* The italics are due to the present author. The reason for this will ap 


pear subsequently 


ferential equation exists is guaranteed by Theorem 3 refereng 
. . . . , eC 
$. It can be exhibited by direct computation 
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On the Boundary-Layer Equations in 
Hypersonic Flow 


S. F. Shen 
Assistant Professor, Department of Aeronautica! Engineering 
University of Maryland, College Park, Md 


April 9, 1952 


HE PROBLEM OF VISCOSITY EFFECTS in hypersonic flow has 
fy pandionns received attention. In 1949, the author made an 
attempt to treat the interaction between shock wave and vis 
cosity near the vertex of an insulated wedge,' using essentially the 
Pohlhausen approach and assuming the whole régime between 
the shock and the body surface describable by the boundary-layer 
equations. The flat plate was solved as a special case. Lees 
and Probstein? solved the complementary problem of the bound 
ary-layer flow over a flat plate at very large distances from the 
leading edge, where interaction may be neglected. The com 
plete Navier-Stokes equations are used. The detail distribu 
tions are obtained by expanding the solutions in ascending powers 
of the parameter j/* V/ Rex, where J/ is the free-stream Mach 
Number and Rex is the Reynolds Number based on distance 4 
from the leading-edge and free-stream conditions. Because of 
this work, it seems that the validity of the boundary-layer equa 
tions in hypersonic flow could be questioned. The present note 
is intended to clarify this situation and to bring out the signifi 
cance of the parameter 17°*/VY Rex. 

It may first be noted that the two problems treated have in 
common the fact that the pressure gradient is essentially due to 
a “‘self-induced”’ effect * depending on the boundary-layer thick 
ness. In the “interaction problem,’’ the pressure gradient is 
given by the oblique shock relations. In the ‘‘boundary-layer 
problem,’’ Lees and Probstein assumed the pressure gradient as 
due to an effective Prandtl-Meyer expansion around the dis- 
placement thickness. Reasonably simple expressions are there- 
fore available for either case. 

However, Tsient pointed out that details of the interaction of 
the flow behind a shock wave and the boundary layer at very 
large Mach Numbers is a complicated flow problem. As indi 
cated in Fig. 1, the shock wave is continuously being weakened 
as the distance from the leading edge increases. The pressurt 
gradient along the wall is thus ‘‘favorable’’ for the boundary 


* This term is used in a letter to the author from R. F. Probstein 
T This paragraph is quoted from a private communication to the author 
by Prof. H.S. Tsien 
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layer 
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SOIT TT 


layer, and there is a region or a layer immediately ‘‘inside’’ of 
the shock where the viscous forces are small but where, never- 
theless, the flow is rotational due to entropy variations generated 
Therefore, beyond the true boundary 


by the curved shock 
In this layer the gas is expand 


laver there is an entropy layer. 
ing However, strictly speaking, the flow in the entropy layer 
cannot be considered as a Prandtl-Meyer flow as is done by Lees 
and Probstein, since Prandtl-Meyer flow is irrotational. In par- 
ticular, it would be incorrect to estimate the pressure gradient 
dp/dy at the edge of the boundary layer as M(0p/0x) 

Navier-Stokes equations, for two-dimensional 


Consider the 
motion, 
ri) re) Op a) Ou Oo Ou 
pu? + puv = t Me t Me 1 
; Oy Ox Ox Ox oy oy 


Ox 


| 


(1) 
where conventional notation has been adopted. To determine 
the proper form of the boundary-layer equation in hypersonic 
flow, one notes 
~ 0 [U"(6/L)] 
~ 0 (6) 


u~ OCU"), v 
x ~ OL), y 
w™ O(n’ M*) 


p~ Op 


with superscript 0 indicating the free-stream quantities; 6, the 


boundary-layer thickness; L, a length scale; and ., the free 
stream Mach Number. The powers \/? and M®* are introduced 
for cases where heat transfer is zero or small at the wall, so that 
the temperature in the boundary layer is of the order of the stag 


nation temperature: 


T’ and Ty being the free-stream and the stagnation temperatures, 
respectively Exponent s is from the viscosity law, u4« 7°. The 
pressure across the boundary layer is implied to remain of the 
same order of magnitude 

Following the usual method of discussion, in boundary layer 


Oo oO Ou 
pu? im ~ (1) 
Ox oy oy 


Substituting the respective orders of magnitude, one gets 


there must be 


6/L ~ 0 (M*t!/¥V R) (2) 
where R = LU°Lp/u’, the Reynolds Number based on free-stream 
The numerator 1/°*', of course, 


conditions and length scale L 
The ‘‘higher 


is due to the hypersonic assumption that 1 > 1 
order” viscosity term (0/Ox) [u(O0u/Ox)|] is, as usual, of the order 


Oo ( 2 re) *) / (:) ss 
ax \" ax oy (. Ox iis Sa sed 
One may next estimate the order of magnitude of the ‘“‘self 
induced” pressure gradient in the case of the “interaction prob 
lem” mentioned above. With simplifications on the oblique 
shock conditions based on the hypersonic assumption, one has! 


FORUM 501 
Op dé d?6 6° 
~ O[ poe ~ 0O( Pur 
Ox dx dx? L3 
hence, 
Op (3) & 
pu2 ~ OT AM? (4a) 
Oox/ O L* 
or by Eq. (2), 
Op 
pur ~ 0 (M R) (4b) 
Ox/ OX 
The following nondimensional variables may now be intro- 
duced: 
ga z'L, y= y'6 
u= Ux’, vy = v' U"(6/L) (5) 
p = p'p M-?, w= pp M 
pb = p'p’|M(6/L)P 


Eq. (1) then takes the form: 


o ne o rpeee o , ou’ 
(pu 2) + _(puv) = (Lew 7 
oy oy oy 


Ox 
2 Op’ u+2 Oo ( , ou’ 
sa Lu -I+... (6) 
Ox R ox Ox 


It is at once 


M?* 
VR 


In Eq. (6), all primed quantities are of order unity. 


that (a) the ‘‘self-induced”’ pressure gradient cannot be 


neglected if 


clear 


|.M(6/L)|? ~ (AMP? VR)? ~ O11) 


(b) the “higher order’ viscosity term, as compared with the 
pressure gradient term, is of the order 


re) ( =) /Op 
M 
Ox Ox// Ox 


It can further be shown from the second momentum equation 


that the pressure gradient across the boundary layer is of the 


~ 0 (M~?) (7) 


order 
(8) 


op’/oy’ ~ 0(M 
Change of p’ across the layer is of the same order as Op’/dy’. 


The hypersonic flow may be defined by 


K ~ M(6/L) ~ O(1) (9) 
Then, both (0/O0x) [u(Ou/dx)| and Op/dy are negligible and the 
ordinary form of the Prandtl boundary-layer equations is still 
valid. 

By taking s = 1, Lees and Probstein arrived at the parameter 
M3/¥V/ Rex, which indeed may be regarded as the same hypersonic 
parameter A defined by Eq. (9). Its significance is no more and 
no less than that of A—namely, defining the hypersonic flow 
On the “boundary-layer” problem, if one took the gradient 0p/0x 
as due to a Prandtl-Meyer flow as in reference 2, similar conclu 
sions regarding the orders of magnitude may be derived. As a 
result, the validity of the boundary-layer equations is not tied in 
with the magnitude of the hypersonic parameter, as long as both 


6/L<land M > 1 hold. 
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Pressure Distribution on an Airfoil in 
Nonuniform Motion 


Robert H. Scanlan* 


Ecole Nationale Supérieure de |’Aéronautique, Paris 
April 10, 1952 


iy A RECENT NOTE in this Forum, Neumark! points out that he 


has ‘“‘succeeded in deriving the long anticipated formula’”’ for 
the pressure distribution on a two-dimensional airfoil in arbi- 
trary motion. 

It may be mentioned that, for the purely sinusoidal case, the 
following more explicit formula already exists in an elementary 
text? and has been derived after the original notions of Theo 
dorsen* following a certain rational approach suggested and 
used by Garrick.‘ 


dp w(x) l x? 
(x, %1,t) = —2pv = 
dx, n(x x1) l x? 
pb w(x), t) (x —)? + [V1 x? Vi x)? ]? 
log 1 
7 ol (x —-x)? + [V1 vr+Vi x? |? 
on 
“pv Ll +2 . F . 
y, w(x, ) 1C(k) +x [1 C(k)]} (1) 
Vil — x? l-—™ 
where x, x; are coordinates on the airfoil between 1 and 1, 


w(x, t) is the value at x; of a vertical velocity “jump” of width 
dx, of the airfoil, v is airfoil forward velocity, and k = bw/v is 
reduced frequency. 

The function p(x, ¢) is given, for a given w-distribution, 


Os v Oz 


w(x; t) = t 
Ox, b ol 


[s(x,) being the vertical deflection distribution of the airfoil] by 
integration over the airfoil with respect to x. 

Postel and Leppert® have explicitly specified the pressure dis 
tribution for the sinusoidal case examined initially by Theodor- 
sen. 

Now it is well known (as has been extensively discussed in the 
March, 1952, Readers’ Forum, for example) that the Theodorsen 
C(k) can be “generalized” for arbitrary specified motions of the 
airfoil by replacing it by a superposition of indicial (Wagner 
function) effects involving a factor consisting of the time deriv- 
ative of airfoil vertical velocity. Specifically, Garrick presents 
the following terms, in the case of completely arbitrary motion, 
to replace only those terms of Eq. (1) having the factor C(k): 

i Py) 
n w(x, O)Ri(s) +n ki(s a) w(x, 0) da (2) 
0 Os 


where 
ks = wt 
2pv l x l+x 
n =n(x,%) = 
T l+vx l x 
and k,(s) is Wagner’s indicial lift function [k)(0) = '/s, ki( ©) = 
1| 


Hence Eq. (1) has considerable generality, since it can be ap 
plied in this manner to unsteady effects by generalization of its 
C(k) terms only, other terms remaining valid as they stand. In 
fact, this generality is effectively superior to that of Neumark’s 
Eq. (2), which contains the important wake contribution as the 
factor, 

¥y(&, t) dé 
ey Ve l 
which appears in both his lift and pressure distribution formulas 


* National Research Fellow, Ecole Nationale Supérieure de 1’ Aéronau 
tique, Paris, France (Associate Professor of Aeronautical Engineering 
Rensselaer Polytechnic Institute, Troy, N.Y., on leave 


This factor, the heart of the problem, requires concrete knowledg 
of the wake vorticity y so that its apparently simple form really 
requires further specification of the relation of the wake to the 
type of motion to be studied. This specification is indeed mad 
» CtC., more 
complicated cases of nonuniform motion being left undone The 


by Neumark for the simple cases of sinusoidal flapping 


results presented in Eqs. (1) and (2) make use of the prior Wagner 
study of the relation of wake vorticity to airfoil motion, ang 


agar 
be solved as is the case with the formulation (3). The simph 


hence, for more complicated motions, this problem need not 


sinusoidal cases worked out by Neumark are, of course, implicit 
in Eq. (1). 

Finally, a ‘historical’? remark is in order on such formulas 5 
Eq. (1). Garrick‘ exhibits an equivalent form and comments as 
follows: 


“The result (p(x, t)] . . . was also explicitly arrived at by 
Kiissner and Schwarz,® by Schwarz,’ and by Séhngen® 

A listing here of two additional and equivalent forms in which 
the result appears in the literature will be useful. The first 
of these is due to Cicala....” 


It is finally noteworthy that Theodorsen*® has performed 
many of the integrations involved in applications of Eq. (1) to 
typical wing motion 


REFERENCES 


Neumark, S., Pressure Distribution on an Airfoil in Nonuniform Moti 
Readers’ Forum, Journal of the Aeronautical Sciences, Vol. 19, No. 13, pp 
214-215, March, 1952 
? Scanlan, R. H., and Rosenbaum, R., Aircraft Vibration and Flutter, pp 
393-394; The Macmillan Company, New York, 1951 
Theodorsen, T., General Theory of Aerodynamic Instability and the Mech 
anism of Flutter, N.A.C.A. T.R. No. 496, 1934 
* Garrick, I. E., Nonsteady Wing Characteristics, Section H, Vol. VI 
Applied High Speed Aerodynamics (Aeronautics Publication Program, Prince 
ton University, 1951), to be published, 1952 
Postel, E. E., and Leppert, E. L., Jr., Theoretical Pressure Distributi 
or Thin Airfoil Oscillating in Incompressible Flow, Journal of the Aeronau 
5, No. 8, pp. 485-492, August, 1948 
* Kiissner, H. G., and Schwarz, L., Luftfahrtforschung 17, pp. 337-354 
1940 (N.A.C.A. T.M. No. 991) 
? Schwarz, L., Luftfahrtforschung 17, pp. 379-386, 1940 
* Sohngen, H., Luftfahrtforschung 17, pp. 401-419, 1940 (British R.T.P 
Tr. 1446 


tical Sciences, Vol 


Analysis of the Elastic and Plastic Stability of 
Sandwich Plates by the Method of Split 
Rigidities—II1 


P. P. Bijlaard 
Professor of Structural Engineering, Cornell University, Ithaca, N.Y 
March 19, 1952 


IT A RECENT PAPER, Seide! recalculated the compressive buck 
ling stress of sandwich plates with clamped unloaded edges 
He finds buckling stresses that, in a certain region, are slightly 
lower than those obtained by the method of split rigidities. It is 
the first case, to the writer’s knowledge, where the method of split 
rigidities overestimates the buckling stress. This requires an 
explanation 
The formula 


Pep = Po + (P,—! + Pe-!) ( 


was derived in reference 2 in a simple way, with the assumption 
that the deflection surfaces w; and w» for Cases 1 (assuming the 
core infinitely rigid against shear) and 2 (assuming only shear de 
formation of the core) had the same shape. However, the same 
formula is found if, in the combined case, the deflection surfaces 


for Cases 1 and 2 differ in shape but have the same shapes as @ 
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nd w2 for the separate individual Cases 1 and 2, respectively. 
This may be proved by comparing the work done by restraining 
und deflecting forces instead of comparing these forces them- 


selves, aS Was suggested in reference 3. This method has since 
been applied successfully to a number of complicated buckling 
problems, of which solutions by other methods seemed to be im 
possible 
As explained in reference 2, the deflecting forces in Case 1 are 
D, = CiPiws With an arbitrary deflection surface C; will be a 
function of the ordinates x and y, which differs from point to 


point and which will be denoted as a, so that D; = a? ,w, per unit 
plate surface. Hence, also, the restraining force in Case | is 


R, = a Pw, (2) 


In the same way the restraining forces R», per unit plate surface, 


generated by the deflection w., which differs in shape from w, are 


R, = BPww> (3) 


where 8 is another function of x and y 

Starting from Case 1, with deflection w,, and using the same 
reasoning as given in reference 2, the restraining forces do not 
increase by the deflection w2, so that for the combined case, with 
+ we, the restraining forces remain /; from 


= W) 


deflection w 
Eq. (2). Analogous to Eqs. (2) and (3) the deflecting forces for 
the combined case are 

Bw ) (4) 
where P?, is the reduced load, as defined in reference 2, which 
neglects the proper flexural rigidity of the faces. With the total 
deflection w = w,; + we, the total negative work done by the re 
straining forces should be equal to the total work done by the de 
flecting forces, so that 

1/2 LS R\(w, + w.) dx dy = 


or, from Eqs. (2) and (4), 


P ime aw, (w, + we) dx dy = 


(1/2) ff D(w, 


+ w»)dx dy 


Pp, SS (cow + BWe)X 


(w; + we) dx dy (5) 


Similarly, analogous to the reasoning in reference 2 and starting 
now from Case 2, the restraining forces for the total deflection re 


main Ry from Eq. (3), so that one obtains 
1/2) [ie i RAw, + w.) dx dy = (1/2) , at 4 D(w, + we) dx dy 
or 
P, ws Bw(w, + w.) dx dy = 
P, ym (aw, + Bw) (wi + we) dx dy (6) 
Multiplication of Eq. (5) by P2 and of Eq. (6) by ?; and addition 
of the results yields 
P,P se Fs (aw, + Bw) (wi + we) dx dy = (P; 4+ 


P, SS (cows 1 Swe) (w, 


P2)X 
+ wy) dx dy 


or 


) 


(P,~' + P.-')* ( 


As in reference 2, addition of the critical thrust Py) of the faces 
leads to the critical load from Eq. (1) according to the method of 


P, = P\P2/(P; + P2) = 


split rigidities 

Hence, Eq. (1) applies if, in the combined case, the deflection 
surfaces w; and w» have the same shape as they have in the com 
ponent Cases 1 and 2, respectively. This is a reasonable assump 
tion and explains why, in cases where w; and we differ in shape, 
relatively accurate results are obtained. 

In reality, the above assumption will not be completely satis 
fied, so that actually in the combined case component deflections 
1’ and wy' occur which differ in shape from w, and w» for the 


w 
individual Cases 1 and 2 Since wz, and we are the optimal de 
flection surfaces, which lead to the individual minimum critical 


loads P; and Ps, the critical loads P;' and P2', belonging to the 
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deflection surfaces w,' and w,', will be higher than ?; and Pz, 
respectively, provided that w,'! and we! have to satisfy the same 
boundary conditions as the combined surface w. Hence, using 


the same reasoning as above, the actual critical load is 
(Perla = Po! + ((Pi')7! + (P32!) ]7 (1’) 


so that, in general, if the deflection surfaces w; and we» have dif- 
ferent shapes, the actual critical load (P.,)q is larger than P,, from 
Eq. (1) as computed by the method of split rigidities. Up to now 
this proved indeed to be true for all cases where comparisons 
with exact calculations could be made. 

However, as mentioned in reference 2, it should be possible 
that the deformations in Cases 1 and 2 occur independently of 
each other. This condition is also satisfied for the present case, 
but it may be shown by sketching the actual deformations that 
the actual Cases 1! and 2! are not quite independent of each other 
Cases 1 and 


2 have to satisfy the same boundary and continuity conditions 


This may also be expressed in the following way. 


as the combined case, but Cases 1' and 2! have to satisfy them 


together only but not separately. This causes a relaxation of 
restraints, so that, although the difference in shape of the deflec- 
tion surfaces w,! and w.!' from w,; and wy tends to make P;' and 
P,' larger than P; and Ps, the former may still be smaller than 
P; and P», so that (P.,)q¢ from Eq. (1’) may be somewhat smaller 
than P,, from Eq. (1). 


to disprove this possibility, from reference 1 this is apparently 


Although Seide’s earlier results seemed 


the case for sandwich plates with clamped unloaded edges. 

This same relaxation of restraints, however, limits the appli 
cability of the formulas for sandwich plates with anistropic core, 
mentioned in the Appendix of reference 4, to ratios ¢ = Gey /Gey 


between about 0.3 and 3. 
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Once More—Single Degree of Freedom Flutter 
of an Aileron 


K. P. Abichandani and R. M. Rosenberg 
Canadair Limited, Montreal, Canada, and Boeing Airplane 
Company, Seattle, Wash., Respectively 
April 16, 1952 
N A CONTRIBUTION TO THIS FORUM, Runyan, Cunningham, 
I and Watkins! discuss a paper (by us) on aileron flutter in the 
single degree of freedom of aileron rotation about the hinge.* 
In our paper, the conclusion was reached that ailerons capable of 
1 with 


be overbalanced It is this conclusion 


this flutter must 
which Runyan, et al., take issue because it is at variance with a 
a result 


result found by these authors* and reported merely as 


and without the theoretical development which, of course, pre 
ceded their finding 

Runyan, et al., state! that we ‘‘concluded that the flutter speed 
the 34 per 


of an aileron hinged ahead of (approximately cent 


aileron chord is imaginary. This conclusion was based, how 


ever, on the (our) failure to realize that the aileron divergence 


speed, rather than the flutter speed, is imaginary for these hinge 
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conditions.’’ Furthermore, they say that with this (i.e., their) 
interpretation, our conclusion ‘‘that the ratio of flutter frequency 
to aileron natural frequency is less than unity should be reversed.”’ 
We shall discuss these points in order. 

To set matters straight with regard to the astonishing con- 
clusions that Runyan, et al., deduced from our paper, we should 
like to point out that our paper is divided into several portions, 
one of which is entitled ‘‘Aileron Divergence.’’ It is in this, 
and only in this, section that we refer to the 34 per cent aileron 
chord. In fact we say that “it has been found that an aileron 
with approximately 34 per cent aerodynamic balance will diverge 
at all speeds.’’ Therefore, we cannot agree that we failed to 
realize that it is the aileron divergence speed that becomes im 
aginary for the hinge location in question (foreward of the 34 
It seems to us that we actually 


per cent aileron chord point). 
found and emphasized what we are said to have failed to realize. 
However, to point out and rectify the fact that Runyan, et al., 
misunderstood and then misquoted our paper does not resolve a 
real discrepancy between the results found by Runyan, Cunning 
ham, and Watkins’ and those found by ourselves. 
quite true that we found large amounts of aerodynamic aileron 


For it is 


balance to be required for flutter, while Runyan, et al., found, on 
the contrary, that ailerons hinged at the leading edge can flutter 
To quote from our paper,” we concluded that ‘‘an aileron capable 
of fluttering (in the sense of this paper) is one which is over 
Our reason for this conclusion is, however, distinct 
If an aileron 


balanced.”’ 
from considerations regarding aileron divergence. 
is to flutter, the system must satisfy simultaneously two condi 
tions [Eqs. (5) of our paper] and no system not satisfying both 
can flutter (in the narrow sense in which the term is used here). 
In particular, the fluttering system must satisfy the condition 
that the imaginary part of the oscillatory aerodynamic torque 
about the hinge must vanish. In our paper, this condition is 
satisfied in Eq. (8), and the results are plotted in Fig. 1. The 
numerical coefficients used in our computations were taken from 
the tables of Smilg and Wasserman ;‘ as a result of these computa 
tions, we find that ailerons hinged at the leading edge do not flut- 
ter. In fact, the equation preceding Eq. (8) in our paper shows 
that, when the aileron aerodynamic balance is zero, ailerons can 
flutter only if the ratio of the imaginary parts of 7, and P; (in 
the nomenclature of reference 4) or 6;/b, (in our nomenclature ) 
vanishes. For v/bw > 0, the tables of Smilg and Wasserman‘ do 
not show any values of 7g and P; which satisfy this criterion. 
An examination of our Fig. 1 will show that, actually, we find flut- 
ter only if the aileron hinge is aft of approximately the 48 per cent 
point of the aileron chord. It is for this reason that we concluded 
that “it is not expected that this type of flutter will be observed in 
flight because it has been shown that only ailerons which are aero- 


dynamically overbalanced can experience it.’”” We have shown? 
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that the divergence speed is always real if the aileron hinge is aft 
of the 34 per cent aileron chord point (approximately) and the 
flutter speed is never real unless the aileron hinge is aft of the 4g 
per cent aileron chord point (approximately). Therefore, the Tas 
tio of flutter to divergence velocities vs/vg cannot be imaginary fog 


real flutter speeds. But this ratio is [Eq. (18), reference 2] 


(v¢/vg)? = [1 — (w/weg)?]/n 


Runyan, Cunningham, and Watkins’ concluded that vy/vg could 
be imaginary because of imaginary vg and real vy. This conclusiog 
is based, however, on the failure of these authors to realize that 
our theoretical development? does not admit this possibility, 
The real disagreement between the results of these authors ang 
ours? remains with considerations of the flutter velocity alone; 
questions of the divergence speed do not enter 

With regard to the point raised by Runyan, Cunningham, and 
Watkins! concerning the ratio w/wg of flutter frequency to aileron 
natural frequency, we suspect that substantial agreement exists, 
provided some slight changes are made in the statements con. 
tained in the March issue of this forum.' We have shown that 
the frequency ratio in question must be less than unity in the 
case of aerodynamically overbalanced ailerons. However, if 
ailerons that are not overbalanced can flutter (vg imaginary, vy 
a frequency greater than 
deduced at 


real!), then they will indeed flutter at 
the aileron natural frequency. This 
once from the above equation, and it is quite in accord with our 
However, this observation does not require the reversal 


fact can be 


findings. 
of any of our conclusions. 

Since the appearance of our paper, a flutter model has been 
constructed in accordance with our paper and tested. The re- 
sults of these tests show gratifying agreement between theory and 
experiment with regard to the flutter speed and the flutter fre- 
In particular, the experiments show that the flutter 


quency. 
However, 


frequency is less than the aileron natural frequency. 
we hasten to point out that these tests do not invalidate (nor do 
they even concern) the results of Runyan, Cunningham, and 
Watkins,’ for no ailerons hinged at the leading edge were tested. 
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